代码

// DUCK_PUBLIC_SURROGATE_BENCHMARK
// Analytic public matrices, independent of production data and all generators.
#pragma GCC optimize("O3")
#pragma GCC target("avx2,fma")
#include <immintrin.h>
#include <stdint.h>
#include <stddef.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
static constexpr int N=8192, ROWS=384, COLS=8160, NC=96, KC=256, MC=48;
alignas(4096) static float bp[NC*KC];
static uint64_t tick(){uint32_t lo,hi;asm volatile("lfence;rdtsc;lfence":"=a"(lo),"=d"(hi)::"memory");return ((uint64_t)hi<<32)|lo;}
namespace m6n16u1 {
constexpr int MR=6,NR=16;
static inline void micro(const float*a,const float*b,float*c,int stride){
 __m256 c00=_mm256_loadu_ps(c+0*stride+0);
 __m256 c01=_mm256_loadu_ps(c+0*stride+8);
 __m256 c10=_mm256_loadu_ps(c+1*stride+0);
 __m256 c11=_mm256_loadu_ps(c+1*stride+8);
 __m256 c20=_mm256_loadu_ps(c+2*stride+0);
 __m256 c21=_mm256_loadu_ps(c+2*stride+8);
 __m256 c30=_mm256_loadu_ps(c+3*stride+0);
 __m256 c31=_mm256_loadu_ps(c+3*stride+8);
 __m256 c40=_mm256_loadu_ps(c+4*stride+0);
 __m256 c41=_mm256_loadu_ps(c+4*stride+8);
 __m256 c50=_mm256_loadu_ps(c+5*stride+0);
 __m256 c51=_mm256_loadu_ps(c+5*stride+8);
#pragma GCC unroll 1
 for(int k=0;k<KC;++k){
 __m256 b0=_mm256_load_ps(b+0);
 __m256 b1=_mm256_load_ps(b+8);
 { __m256 a0=_mm256_broadcast_ss(a+0);c00=_mm256_fmadd_ps(a0,b0,c00);c01=_mm256_fmadd_ps(a0,b1,c01);}
 { __m256 a1=_mm256_broadcast_ss(a+1);c10=_mm256_fmadd_ps(a1,b0,c10);c11=_mm256_fmadd_ps(a1,b1,c11);}
 { __m256 a2=_mm256_broadcast_ss(a+2);c20=_mm256_fmadd_ps(a2,b0,c20);c21=_mm256_fmadd_ps(a2,b1,c21);}
 { __m256 a3=_mm256_broadcast_ss(a+3);c30=_mm256_fmadd_ps(a3,b0,c30);c31=_mm256_fmadd_ps(a3,b1,c31);}
 { __m256 a4=_mm256_broadcast_ss(a+4);c40=_mm256_fmadd_ps(a4,b0,c40);c41=_mm256_fmadd_ps(a4,b1,c41);}
 { __m256 a5=_mm256_broadcast_ss(a+5);c50=_mm256_fmadd_ps(a5,b0,c50);c51=_mm256_fmadd_ps(a5,b1,c51);}
 a+=MR;b+=NR;
 }
 _mm256_storeu_ps(c+0*stride+0,c00);
 _mm256_storeu_ps(c+0*stride+8,c01);
 _mm256_storeu_ps(c+1*stride+0,c10);
 _mm256_storeu_ps(c+1*stride+8,c11);
 _mm256_storeu_ps(c+2*stride+0,c20);
 _mm256_storeu_ps(c+2*stride+8,c21);
 _mm256_storeu_ps(c+3*stride+0,c30);
 _mm256_storeu_ps(c+3*stride+8,c31);
 _mm256_storeu_ps(c+4*stride+0,c40);
 _mm256_storeu_ps(c+4*stride+8,c41);
 _mm256_storeu_ps(c+5*stride+0,c50);
 _mm256_storeu_ps(c+5*stride+8,c51);
}
static void run(const float*A,const float*B,float*C){
 memset(C,0,(size_t)ROWS*N*4);
 float*ap=nullptr; if(posix_memalign((void**)&ap,4096,(size_t)ROWS*N*4))abort();
 constexpr int groups=ROWS/MR;
 for(int pc=0;pc<N;pc+=KC)for(int g=0;g<groups;++g)for(int k=0;k<KC;++k)for(int i=0;i<MR;++i)
  ap[(((size_t)(pc/KC)*groups+g)*KC+k)*MR+i]=A[(size_t)(g*MR+i)*N+pc+k];
 uint64_t begin=tick(),packing=0;
 for(int jc=0;jc<COLS;jc+=NC)for(int pc=0;pc<N;pc+=KC){
  uint64_t ps=tick();
  for(int q=0;q<NC;q+=NR)for(int k=0;k<KC;++k)for(int j=0;j<NR;j+=8)
   _mm256_store_ps(bp+((q/NR*KC+k)*NR)+j,_mm256_loadu_ps(B+(size_t)(pc+k)*N+jc+q+j));
  packing+=tick()-ps;
  for(int ic=0;ic<ROWS;ic+=MC)for(int q=0;q<NC;q+=NR)for(int r=0;r<MC;r+=MR)
   micro(ap+((size_t)(pc/KC)*groups+(ic+r)/MR)*KC*MR,bp+q*KC,C+(size_t)(ic+r)*N+jc+q,N);
 }
 uint64_t total=tick()-begin;
 for(int t=0;t<4;++t){int i=t*91,j=t*1999;double expected=0;
  for(int k=0;k<N;++k)expected+=(double)A[(size_t)i*N+k]*B[(size_t)k*N+j];
  if(fabs(expected-C[(size_t)i*N+j])>0.1)abort();
 }
 printf("m6n16u1 total_ticks=%llu pack_ticks=%llu core_ticks=%llu\n",(unsigned long long)total,(unsigned long long)packing,(unsigned long long)(total-packing));
 free(ap);
}
}
namespace m6n16u2 {
constexpr int MR=6,NR=16;
static inline void micro(const float*a,const float*b,float*c,int stride){
 __m256 c00=_mm256_loadu_ps(c+0*stride+0);
 __m256 c01=_mm256_loadu_ps(c+0*stride+8);
 __m256 c10=_mm256_loadu_ps(c+1*stride+0);
 __m256 c11=_mm256_loadu_ps(c+1*stride+8);
 __m256 c20=_mm256_loadu_ps(c+2*stride+0);
 __m256 c21=_mm256_loadu_ps(c+2*stride+8);
 __m256 c30=_mm256_loadu_ps(c+3*stride+0);
 __m256 c31=_mm256_loadu_ps(c+3*stride+8);
 __m256 c40=_mm256_loadu_ps(c+4*stride+0);
 __m256 c41=_mm256_loadu_ps(c+4*stride+8);
 __m256 c50=_mm256_loadu_ps(c+5*stride+0);
 __m256 c51=_mm256_loadu_ps(c+5*stride+8);
#pragma GCC unroll 2
 for(int k=0;k<KC;++k){
 __m256 b0=_mm256_load_ps(b+0);
 __m256 b1=_mm256_load_ps(b+8);
 { __m256 a0=_mm256_broadcast_ss(a+0);c00=_mm256_fmadd_ps(a0,b0,c00);c01=_mm256_fmadd_ps(a0,b1,c01);}
 { __m256 a1=_mm256_broadcast_ss(a+1);c10=_mm256_fmadd_ps(a1,b0,c10);c11=_mm256_fmadd_ps(a1,b1,c11);}
 { __m256 a2=_mm256_broadcast_ss(a+2);c20=_mm256_fmadd_ps(a2,b0,c20);c21=_mm256_fmadd_ps(a2,b1,c21);}
 { __m256 a3=_mm256_broadcast_ss(a+3);c30=_mm256_fmadd_ps(a3,b0,c30);c31=_mm256_fmadd_ps(a3,b1,c31);}
 { __m256 a4=_mm256_broadcast_ss(a+4);c40=_mm256_fmadd_ps(a4,b0,c40);c41=_mm256_fmadd_ps(a4,b1,c41);}
 { __m256 a5=_mm256_broadcast_ss(a+5);c50=_mm256_fmadd_ps(a5,b0,c50);c51=_mm256_fmadd_ps(a5,b1,c51);}
 a+=MR;b+=NR;
 }
 _mm256_storeu_ps(c+0*stride+0,c00);
 _mm256_storeu_ps(c+0*stride+8,c01);
 _mm256_storeu_ps(c+1*stride+0,c10);
 _mm256_storeu_ps(c+1*stride+8,c11);
 _mm256_storeu_ps(c+2*stride+0,c20);
 _mm256_storeu_ps(c+2*stride+8,c21);
 _mm256_storeu_ps(c+3*stride+0,c30);
 _mm256_storeu_ps(c+3*stride+8,c31);
 _mm256_storeu_ps(c+4*stride+0,c40);
 _mm256_storeu_ps(c+4*stride+8,c41);
 _mm256_storeu_ps(c+5*stride+0,c50);
 _mm256_storeu_ps(c+5*stride+8,c51);
}
static void run(const float*A,const float*B,float*C){
 memset(C,0,(size_t)ROWS*N*4);
 float*ap=nullptr; if(posix_memalign((void**)&ap,4096,(size_t)ROWS*N*4))abort();
 constexpr int groups=ROWS/MR;
 for(int pc=0;pc<N;pc+=KC)for(int g=0;g<groups;++g)for(int k=0;k<KC;++k)for(int i=0;i<MR;++i)
  ap[(((size_t)(pc/KC)*groups+g)*KC+k)*MR+i]=A[(size_t)(g*MR+i)*N+pc+k];
 uint64_t begin=tick(),packing=0;
 for(int jc=0;jc<COLS;jc+=NC)for(int pc=0;pc<N;pc+=KC){
  uint64_t ps=tick();
  for(int q=0;q<NC;q+=NR)for(int k=0;k<KC;++k)for(int j=0;j<NR;j+=8)
   _mm256_store_ps(bp+((q/NR*KC+k)*NR)+j,_mm256_loadu_ps(B+(size_t)(pc+k)*N+jc+q+j));
  packing+=tick()-ps;
  for(int ic=0;ic<ROWS;ic+=MC)for(int q=0;q<NC;q+=NR)for(int r=0;r<MC;r+=MR)
   micro(ap+((size_t)(pc/KC)*groups+(ic+r)/MR)*KC*MR,bp+q*KC,C+(size_t)(ic+r)*N+jc+q,N);
 }
 uint64_t total=tick()-begin;
 for(int t=0;t<4;++t){int i=t*91,j=t*1999;double expected=0;
  for(int k=0;k<N;++k)expected+=(double)A[(size_t)i*N+k]*B[(size_t)k*N+j];
  if(fabs(expected-C[(size_t)i*N+j])>0.1)abort();
 }
 printf("m6n16u2 total_ticks=%llu pack_ticks=%llu core_ticks=%llu\n",(unsigned long long)total,(unsigned long long)packing,(unsigned long long)(total-packing));
 free(ap);
}
}
namespace m6n16u4 {
constexpr int MR=6,NR=16;
static inline void micro(const float*a,const float*b,float*c,int stride){
 __m256 c00=_mm256_loadu_ps(c+0*stride+0);
 __m256 c01=_mm256_loadu_ps(c+0*stride+8);
 __m256 c10=_mm256_loadu_ps(c+1*stride+0);
 __m256 c11=_mm256_loadu_ps(c+1*stride+8);
 __m256 c20=_mm256_loadu_ps(c+2*stride+0);
 __m256 c21=_mm256_loadu_ps(c+2*stride+8);
 __m256 c30=_mm256_loadu_ps(c+3*stride+0);
 __m256 c31=_mm256_loadu_ps(c+3*stride+8);
 __m256 c40=_mm256_loadu_ps(c+4*stride+0);
 __m256 c41=_mm256_loadu_ps(c+4*stride+8);
 __m256 c50=_mm256_loadu_ps(c+5*stride+0);
 __m256 c51=_mm256_loadu_ps(c+5*stride+8);
#pragma GCC unroll 4
 for(int k=0;k<KC;++k){
 __m256 b0=_mm256_load_ps(b+0);
 __m256 b1=_mm256_load_ps(b+8);
 { __m256 a0=_mm256_broadcast_ss(a+0);c00=_mm256_fmadd_ps(a0,b0,c00);c01=_mm256_fmadd_ps(a0,b1,c01);}
 { __m256 a1=_mm256_broadcast_ss(a+1);c10=_mm256_fmadd_ps(a1,b0,c10);c11=_mm256_fmadd_ps(a1,b1,c11);}
 { __m256 a2=_mm256_broadcast_ss(a+2);c20=_mm256_fmadd_ps(a2,b0,c20);c21=_mm256_fmadd_ps(a2,b1,c21);}
 { __m256 a3=_mm256_broadcast_ss(a+3);c30=_mm256_fmadd_ps(a3,b0,c30);c31=_mm256_fmadd_ps(a3,b1,c31);}
 { __m256 a4=_mm256_broadcast_ss(a+4);c40=_mm256_fmadd_ps(a4,b0,c40);c41=_mm256_fmadd_ps(a4,b1,c41);}
 { __m256 a5=_mm256_broadcast_ss(a+5);c50=_mm256_fmadd_ps(a5,b0,c50);c51=_mm256_fmadd_ps(a5,b1,c51);}
 a+=MR;b+=NR;
 }
 _mm256_storeu_ps(c+0*stride+0,c00);
 _mm256_storeu_ps(c+0*stride+8,c01);
 _mm256_storeu_ps(c+1*stride+0,c10);
 _mm256_storeu_ps(c+1*stride+8,c11);
 _mm256_storeu_ps(c+2*stride+0,c20);
 _mm256_storeu_ps(c+2*stride+8,c21);
 _mm256_storeu_ps(c+3*stride+0,c30);
 _mm256_storeu_ps(c+3*stride+8,c31);
 _mm256_storeu_ps(c+4*stride+0,c40);
 _mm256_storeu_ps(c+4*stride+8,c41);
 _mm256_storeu_ps(c+5*stride+0,c50);
 _mm256_storeu_ps(c+5*stride+8,c51);
}
static void run(const float*A,const float*B,float*C){
 memset(C,0,(size_t)ROWS*N*4);
 float*ap=nullptr; if(posix_memalign((void**)&ap,4096,(size_t)ROWS*N*4))abort();
 constexpr int groups=ROWS/MR;
 for(int pc=0;pc<N;pc+=KC)for(int g=0;g<groups;++g)for(int k=0;k<KC;++k)for(int i=0;i<MR;++i)
  ap[(((size_t)(pc/KC)*groups+g)*KC+k)*MR+i]=A[(size_t)(g*MR+i)*N+pc+k];
 uint64_t begin=tick(),packing=0;
 for(int jc=0;jc<COLS;jc+=NC)for(int pc=0;pc<N;pc+=KC){
  uint64_t ps=tick();
  for(int q=0;q<NC;q+=NR)for(int k=0;k<KC;++k)for(int j=0;j<NR;j+=8)
   _mm256_store_ps(bp+((q/NR*KC+k)*NR)+j,_mm256_loadu_ps(B+(size_t)(pc+k)*N+jc+q+j));
  packing+=tick()-ps;
  for(int ic=0;ic<ROWS;ic+=MC)for(int q=0;q<NC;q+=NR)for(int r=0;r<MC;r+=MR)
   micro(ap+((size_t)(pc/KC)*groups+(ic+r)/MR)*KC*MR,bp+q*KC,C+(size_t)(ic+r)*N+jc+q,N);
 }
 uint64_t total=tick()-begin;
 for(int t=0;t<4;++t){int i=t*91,j=t*1999;double expected=0;
  for(int k=0;k<N;++k)expected+=(double)A[(size_t)i*N+k]*B[(size_t)k*N+j];
  if(fabs(expected-C[(size_t)i*N+j])>0.1)abort();
 }
 printf("m6n16u4 total_ticks=%llu pack_ticks=%llu core_ticks=%llu\n",(unsigned long long)total,(unsigned long long)packing,(unsigned long long)(total-packing));
 free(ap);
}
}
namespace m6n16u8 {
constexpr int MR=6,NR=16;
static inline void micro(const float*a,const float*b,float*c,int stride){
 __m256 c00=_mm256_loadu_ps(c+0*stride+0);
 __m256 c01=_mm256_loadu_ps(c+0*stride+8);
 __m256 c10=_mm256_loadu_ps(c+1*stride+0);
 __m256 c11=_mm256_loadu_ps(c+1*stride+8);
 __m256 c20=_mm256_loadu_ps(c+2*stride+0);
 __m256 c21=_mm256_loadu_ps(c+2*stride+8);
 __m256 c30=_mm256_loadu_ps(c+3*stride+0);
 __m256 c31=_mm256_loadu_ps(c+3*stride+8);
 __m256 c40=_mm256_loadu_ps(c+4*stride+0);
 __m256 c41=_mm256_loadu_ps(c+4*stride+8);
 __m256 c50=_mm256_loadu_ps(c+5*stride+0);
 __m256 c51=_mm256_loadu_ps(c+5*stride+8);
#pragma GCC unroll 8
 for(int k=0;k<KC;++k){
 __m256 b0=_mm256_load_ps(b+0);
 __m256 b1=_mm256_load_ps(b+8);
 { __m256 a0=_mm256_broadcast_ss(a+0);c00=_mm256_fmadd_ps(a0,b0,c00);c01=_mm256_fmadd_ps(a0,b1,c01);}
 { __m256 a1=_mm256_broadcast_ss(a+1);c10=_mm256_fmadd_ps(a1,b0,c10);c11=_mm256_fmadd_ps(a1,b1,c11);}
 { __m256 a2=_mm256_broadcast_ss(a+2);c20=_mm256_fmadd_ps(a2,b0,c20);c21=_mm256_fmadd_ps(a2,b1,c21);}
 { __m256 a3=_mm256_broadcast_ss(a+3);c30=_mm256_fmadd_ps(a3,b0,c30);c31=_mm256_fmadd_ps(a3,b1,c31);}
 { __m256 a4=_mm256_broadcast_ss(a+4);c40=_mm256_fmadd_ps(a4,b0,c40);c41=_mm256_fmadd_ps(a4,b1,c41);}
 { __m256 a5=_mm256_broadcast_ss(a+5);c50=_mm256_fmadd_ps(a5,b0,c50);c51=_mm256_fmadd_ps(a5,b1,c51);}
 a+=MR;b+=NR;
 }
 _mm256_storeu_ps(c+0*stride+0,c00);
 _mm256_storeu_ps(c+0*stride+8,c01);
 _mm256_storeu_ps(c+1*stride+0,c10);
 _mm256_storeu_ps(c+1*stride+8,c11);
 _mm256_storeu_ps(c+2*stride+0,c20);
 _mm256_storeu_ps(c+2*stride+8,c21);
 _mm256_storeu_ps(c+3*stride+0,c30);
 _mm256_storeu_ps(c+3*stride+8,c31);
 _mm256_storeu_ps(c+4*stride+0,c40);
 _mm256_storeu_ps(c+4*stride+8,c41);
 _mm256_storeu_ps(c+5*stride+0,c50);
 _mm256_storeu_ps(c+5*stride+8,c51);
}
static void run(const float*A,const float*B,float*C){
 memset(C,0,(size_t)ROWS*N*4);
 float*ap=nullptr; if(posix_memalign((void**)&ap,4096,(size_t)ROWS*N*4))abort();
 constexpr int groups=ROWS/MR;
 for(int pc=0;pc<N;pc+=KC)for(int g=0;g<groups;++g)for(int k=0;k<KC;++k)for(int i=0;i<MR;++i)
  ap[(((size_t)(pc/KC)*groups+g)*KC+k)*MR+i]=A[(size_t)(g*MR+i)*N+pc+k];
 uint64_t begin=tick(),packing=0;
 for(int jc=0;jc<COLS;jc+=NC)for(int pc=0;pc<N;pc+=KC){
  uint64_t ps=tick();
  for(int q=0;q<NC;q+=NR)for(int k=0;k<KC;++k)for(int j=0;j<NR;j+=8)
   _mm256_store_ps(bp+((q/NR*KC+k)*NR)+j,_mm256_loadu_ps(B+(size_t)(pc+k)*N+jc+q+j));
  packing+=tick()-ps;
  for(int ic=0;ic<ROWS;ic+=MC)for(int q=0;q<NC;q+=NR)for(int r=0;r<MC;r+=MR)
   micro(ap+((size_t)(pc/KC)*groups+(ic+r)/MR)*KC*MR,bp+q*KC,C+(size_t)(ic+r)*N+jc+q,N);
 }
 uint64_t total=tick()-begin;
 for(int t=0;t<4;++t){int i=t*91,j=t*1999;double expected=0;
  for(int k=0;k<N;++k)expected+=(double)A[(size_t)i*N+k]*B[(size_t)k*N+j];
  if(fabs(expected-C[(size_t)i*N+j])>0.1)abort();
 }
 printf("m6n16u8 total_ticks=%llu pack_ticks=%llu core_ticks=%llu\n",(unsigned long long)total,(unsigned long long)packing,(unsigned long long)(total-packing));
 free(ap);
}
}
namespace m4n24u1 {
constexpr int MR=4,NR=24;
static inline void micro(const float*a,const float*b,float*c,int stride){
 __m256 c00=_mm256_loadu_ps(c+0*stride+0);
 __m256 c01=_mm256_loadu_ps(c+0*stride+8);
 __m256 c02=_mm256_loadu_ps(c+0*stride+16);
 __m256 c10=_mm256_loadu_ps(c+1*stride+0);
 __m256 c11=_mm256_loadu_ps(c+1*stride+8);
 __m256 c12=_mm256_loadu_ps(c+1*stride+16);
 __m256 c20=_mm256_loadu_ps(c+2*stride+0);
 __m256 c21=_mm256_loadu_ps(c+2*stride+8);
 __m256 c22=_mm256_loadu_ps(c+2*stride+16);
 __m256 c30=_mm256_loadu_ps(c+3*stride+0);
 __m256 c31=_mm256_loadu_ps(c+3*stride+8);
 __m256 c32=_mm256_loadu_ps(c+3*stride+16);
#pragma GCC unroll 1
 for(int k=0;k<KC;++k){
 __m256 b0=_mm256_load_ps(b+0);
 __m256 b1=_mm256_load_ps(b+8);
 __m256 b2=_mm256_load_ps(b+16);
 { __m256 a0=_mm256_broadcast_ss(a+0);c00=_mm256_fmadd_ps(a0,b0,c00);c01=_mm256_fmadd_ps(a0,b1,c01);c02=_mm256_fmadd_ps(a0,b2,c02);}
 { __m256 a1=_mm256_broadcast_ss(a+1);c10=_mm256_fmadd_ps(a1,b0,c10);c11=_mm256_fmadd_ps(a1,b1,c11);c12=_mm256_fmadd_ps(a1,b2,c12);}
 { __m256 a2=_mm256_broadcast_ss(a+2);c20=_mm256_fmadd_ps(a2,b0,c20);c21=_mm256_fmadd_ps(a2,b1,c21);c22=_mm256_fmadd_ps(a2,b2,c22);}
 { __m256 a3=_mm256_broadcast_ss(a+3);c30=_mm256_fmadd_ps(a3,b0,c30);c31=_mm256_fmadd_ps(a3,b1,c31);c32=_mm256_fmadd_ps(a3,b2,c32);}
 a+=MR;b+=NR;
 }
 _mm256_storeu_ps(c+0*stride+0,c00);
 _mm256_storeu_ps(c+0*stride+8,c01);
 _mm256_storeu_ps(c+0*stride+16,c02);
 _mm256_storeu_ps(c+1*stride+0,c10);
 _mm256_storeu_ps(c+1*stride+8,c11);
 _mm256_storeu_ps(c+1*stride+16,c12);
 _mm256_storeu_ps(c+2*stride+0,c20);
 _mm256_storeu_ps(c+2*stride+8,c21);
 _mm256_storeu_ps(c+2*stride+16,c22);
 _mm256_storeu_ps(c+3*stride+0,c30);
 _mm256_storeu_ps(c+3*stride+8,c31);
 _mm256_storeu_ps(c+3*stride+16,c32);
}
static void run(const float*A,const float*B,float*C){
 memset(C,0,(size_t)ROWS*N*4);
 float*ap=nullptr; if(posix_memalign((void**)&ap,4096,(size_t)ROWS*N*4))abort();
 constexpr int groups=ROWS/MR;
 for(int pc=0;pc<N;pc+=KC)for(int g=0;g<groups;++g)for(int k=0;k<KC;++k)for(int i=0;i<MR;++i)
  ap[(((size_t)(pc/KC)*groups+g)*KC+k)*MR+i]=A[(size_t)(g*MR+i)*N+pc+k];
 uint64_t begin=tick(),packing=0;
 for(int jc=0;jc<COLS;jc+=NC)for(int pc=0;pc<N;pc+=KC){
  uint64_t ps=tick();
  for(int q=0;q<NC;q+=NR)for(int k=0;k<KC;++k)for(int j=0;j<NR;j+=8)
   _mm256_store_ps(bp+((q/NR*KC+k)*NR)+j,_mm256_loadu_ps(B+(size_t)(pc+k)*N+jc+q+j));
  packing+=tick()-ps;
  for(int ic=0;ic<ROWS;ic+=MC)for(int q=0;q<NC;q+=NR)for(int r=0;r<MC;r+=MR)
   micro(ap+((size_t)(pc/KC)*groups+(ic+r)/MR)*KC*MR,bp+q*KC,C+(size_t)(ic+r)*N+jc+q,N);
 }
 uint64_t total=tick()-begin;
 for(int t=0;t<4;++t){int i=t*91,j=t*1999;double expected=0;
  for(int k=0;k<N;++k)expected+=(double)A[(size_t)i*N+k]*B[(size_t)k*N+j];
  if(fabs(expected-C[(size_t)i*N+j])>0.1)abort();
 }
 printf("m4n24u1 total_ticks=%llu pack_ticks=%llu core_ticks=%llu\n",(unsigned long long)total,(unsigned long long)packing,(unsigned long long)(total-packing));
 free(ap);
}
}
namespace m4n24u2 {
constexpr int MR=4,NR=24;
static inline void micro(const float*a,const float*b,float*c,int stride){
 __m256 c00=_mm256_loadu_ps(c+0*stride+0);
 __m256 c01=_mm256_loadu_ps(c+0*stride+8);
 __m256 c02=_mm256_loadu_ps(c+0*stride+16);
 __m256 c10=_mm256_loadu_ps(c+1*stride+0);
 __m256 c11=_mm256_loadu_ps(c+1*stride+8);
 __m256 c12=_mm256_loadu_ps(c+1*stride+16);
 __m256 c20=_mm256_loadu_ps(c+2*stride+0);
 __m256 c21=_mm256_loadu_ps(c+2*stride+8);
 __m256 c22=_mm256_loadu_ps(c+2*stride+16);
 __m256 c30=_mm256_loadu_ps(c+3*stride+0);
 __m256 c31=_mm256_loadu_ps(c+3*stride+8);
 __m256 c32=_mm256_loadu_ps(c+3*stride+16);
#pragma GCC unroll 2
 for(int k=0;k<KC;++k){
 __m256 b0=_mm256_load_ps(b+0);
 __m256 b1=_mm256_load_ps(b+8);
 __m256 b2=_mm256_load_ps(b+16);
 { __m256 a0=_mm256_broadcast_ss(a+0);c00=_mm256_fmadd_ps(a0,b0,c00);c01=_mm256_fmadd_ps(a0,b1,c01);c02=_mm256_fmadd_ps(a0,b2,c02);}
 { __m256 a1=_mm256_broadcast_ss(a+1);c10=_mm256_fmadd_ps(a1,b0,c10);c11=_mm256_fmadd_ps(a1,b1,c11);c12=_mm256_fmadd_ps(a1,b2,c12);}
 { __m256 a2=_mm256_broadcast_ss(a+2);c20=_mm256_fmadd_ps(a2,b0,c20);c21=_mm256_fmadd_ps(a2,b1,c21);c22=_mm256_fmadd_ps(a2,b2,c22);}
 { __m256 a3=_mm256_broadcast_ss(a+3);c30=_mm256_fmadd_ps(a3,b0,c30);c31=_mm256_fmadd_ps(a3,b1,c31);c32=_mm256_fmadd_ps(a3,b2,c32);}
 a+=MR;b+=NR;
 }
 _mm256_storeu_ps(c+0*stride+0,c00);
 _mm256_storeu_ps(c+0*stride+8,c01);
 _mm256_storeu_ps(c+0*stride+16,c02);
 _mm256_storeu_ps(c+1*stride+0,c10);
 _mm256_storeu_ps(c+1*stride+8,c11);
 _mm256_storeu_ps(c+1*stride+16,c12);
 _mm256_storeu_ps(c+2*stride+0,c20);
 _mm256_storeu_ps(c+2*stride+8,c21);
 _mm256_storeu_ps(c+2*stride+16,c22);
 _mm256_storeu_ps(c+3*stride+0,c30);
 _mm256_storeu_ps(c+3*stride+8,c31);
 _mm256_storeu_ps(c+3*stride+16,c32);
}
static void run(const float*A,const float*B,float*C){
 memset(C,0,(size_t)ROWS*N*4);
 float*ap=nullptr; if(posix_memalign((void**)&ap,4096,(size_t)ROWS*N*4))abort();
 constexpr int groups=ROWS/MR;
 for(int pc=0;pc<N;pc+=KC)for(int g=0;g<groups;++g)for(int k=0;k<KC;++k)for(int i=0;i<MR;++i)
  ap[(((size_t)(pc/KC)*groups+g)*KC+k)*MR+i]=A[(size_t)(g*MR+i)*N+pc+k];
 uint64_t begin=tick(),packing=0;
 for(int jc=0;jc<COLS;jc+=NC)for(int pc=0;pc<N;pc+=KC){
  uint64_t ps=tick();
  for(int q=0;q<NC;q+=NR)for(int k=0;k<KC;++k)for(int j=0;j<NR;j+=8)
   _mm256_store_ps(bp+((q/NR*KC+k)*NR)+j,_mm256_loadu_ps(B+(size_t)(pc+k)*N+jc+q+j));
  packing+=tick()-ps;
  for(int ic=0;ic<ROWS;ic+=MC)for(int q=0;q<NC;q+=NR)for(int r=0;r<MC;r+=MR)
   micro(ap+((size_t)(pc/KC)*groups+(ic+r)/MR)*KC*MR,bp+q*KC,C+(size_t)(ic+r)*N+jc+q,N);
 }
 uint64_t total=tick()-begin;
 for(int t=0;t<4;++t){int i=t*91,j=t*1999;double expected=0;
  for(int k=0;k<N;++k)expected+=(double)A[(size_t)i*N+k]*B[(size_t)k*N+j];
  if(fabs(expected-C[(size_t)i*N+j])>0.1)abort();
 }
 printf("m4n24u2 total_ticks=%llu pack_ticks=%llu core_ticks=%llu\n",(unsigned long long)total,(unsigned long long)packing,(unsigned long long)(total-packing));
 free(ap);
}
}
namespace m4n24u4 {
constexpr int MR=4,NR=24;
static inline void micro(const float*a,const float*b,float*c,int stride){
 __m256 c00=_mm256_loadu_ps(c+0*stride+0);
 __m256 c01=_mm256_loadu_ps(c+0*stride+8);
 __m256 c02=_mm256_loadu_ps(c+0*stride+16);
 __m256 c10=_mm256_loadu_ps(c+1*stride+0);
 __m256 c11=_mm256_loadu_ps(c+1*stride+8);
 __m256 c12=_mm256_loadu_ps(c+1*stride+16);
 __m256 c20=_mm256_loadu_ps(c+2*stride+0);
 __m256 c21=_mm256_loadu_ps(c+2*stride+8);
 __m256 c22=_mm256_loadu_ps(c+2*stride+16);
 __m256 c30=_mm256_loadu_ps(c+3*stride+0);
 __m256 c31=_mm256_loadu_ps(c+3*stride+8);
 __m256 c32=_mm256_loadu_ps(c+3*stride+16);
#pragma GCC unroll 4
 for(int k=0;k<KC;++k){
 __m256 b0=_mm256_load_ps(b+0);
 __m256 b1=_mm256_load_ps(b+8);
 __m256 b2=_mm256_load_ps(b+16);
 { __m256 a0=_mm256_broadcast_ss(a+0);c00=_mm256_fmadd_ps(a0,b0,c00);c01=_mm256_fmadd_ps(a0,b1,c01);c02=_mm256_fmadd_ps(a0,b2,c02);}
 { __m256 a1=_mm256_broadcast_ss(a+1);c10=_mm256_fmadd_ps(a1,b0,c10);c11=_mm256_fmadd_ps(a1,b1,c11);c12=_mm256_fmadd_ps(a1,b2,c12);}
 { __m256 a2=_mm256_broadcast_ss(a+2);c20=_mm256_fmadd_ps(a2,b0,c20);c21=_mm256_fmadd_ps(a2,b1,c21);c22=_mm256_fmadd_ps(a2,b2,c22);}
 { __m256 a3=_mm256_broadcast_ss(a+3);c30=_mm256_fmadd_ps(a3,b0,c30);c31=_mm256_fmadd_ps(a3,b1,c31);c32=_mm256_fmadd_ps(a3,b2,c32);}
 a+=MR;b+=NR;
 }
 _mm256_storeu_ps(c+0*stride+0,c00);
 _mm256_storeu_ps(c+0*stride+8,c01);
 _mm256_storeu_ps(c+0*stride+16,c02);
 _mm256_storeu_ps(c+1*stride+0,c10);
 _mm256_storeu_ps(c+1*stride+8,c11);
 _mm256_storeu_ps(c+1*stride+16,c12);
 _mm256_storeu_ps(c+2*stride+0,c20);
 _mm256_storeu_ps(c+2*stride+8,c21);
 _mm256_storeu_ps(c+2*stride+16,c22);
 _mm256_storeu_ps(c+3*stride+0,c30);
 _mm256_storeu_ps(c+3*stride+8,c31);
 _mm256_storeu_ps(c+3*stride+16,c32);
}
static void run(const float*A,const float*B,float*C){
 memset(C,0,(size_t)ROWS*N*4);
 float*ap=nullptr; if(posix_memalign((void**)&ap,4096,(size_t)ROWS*N*4))abort();
 constexpr int groups=ROWS/MR;
 for(int pc=0;pc<N;pc+=KC)for(int g=0;g<groups;++g)for(int k=0;k<KC;++k)for(int i=0;i<MR;++i)
  ap[(((size_t)(pc/KC)*groups+g)*KC+k)*MR+i]=A[(size_t)(g*MR+i)*N+pc+k];
 uint64_t begin=tick(),packing=0;
 for(int jc=0;jc<COLS;jc+=NC)for(int pc=0;pc<N;pc+=KC){
  uint64_t ps=tick();
  for(int q=0;q<NC;q+=NR)for(int k=0;k<KC;++k)for(int j=0;j<NR;j+=8)
   _mm256_store_ps(bp+((q/NR*KC+k)*NR)+j,_mm256_loadu_ps(B+(size_t)(pc+k)*N+jc+q+j));
  packing+=tick()-ps;
  for(int ic=0;ic<ROWS;ic+=MC)for(int q=0;q<NC;q+=NR)for(int r=0;r<MC;r+=MR)
   micro(ap+((size_t)(pc/KC)*groups+(ic+r)/MR)*KC*MR,bp+q*KC,C+(size_t)(ic+r)*N+jc+q,N);
 }
 uint64_t total=tick()-begin;
 for(int t=0;t<4;++t){int i=t*91,j=t*1999;double expected=0;
  for(int k=0;k<N;++k)expected+=(double)A[(size_t)i*N+k]*B[(size_t)k*N+j];
  if(fabs(expected-C[(size_t)i*N+j])>0.1)abort();
 }
 printf("m4n24u4 total_ticks=%llu pack_ticks=%llu core_ticks=%llu\n",(unsigned long long)total,(unsigned long long)packing,(unsigned long long)(total-packing));
 free(ap);
}
}
namespace m8n8u1 {
constexpr int MR=8,NR=8;
static inline void micro(const float*a,const float*b,float*c,int stride){
 __m256 c00=_mm256_loadu_ps(c+0*stride+0);
 __m256 c10=_mm256_loadu_ps(c+1*stride+0);
 __m256 c20=_mm256_loadu_ps(c+2*stride+0);
 __m256 c30=_mm256_loadu_ps(c+3*stride+0);
 __m256 c40=_mm256_loadu_ps(c+4*stride+0);
 __m256 c50=_mm256_loadu_ps(c+5*stride+0);
 __m256 c60=_mm256_loadu_ps(c+6*stride+0);
 __m256 c70=_mm256_loadu_ps(c+7*stride+0);
#pragma GCC unroll 1
 for(int k=0;k<KC;++k){
 __m256 b0=_mm256_load_ps(b+0);
 { __m256 a0=_mm256_broadcast_ss(a+0);c00=_mm256_fmadd_ps(a0,b0,c00);}
 { __m256 a1=_mm256_broadcast_ss(a+1);c10=_mm256_fmadd_ps(a1,b0,c10);}
 { __m256 a2=_mm256_broadcast_ss(a+2);c20=_mm256_fmadd_ps(a2,b0,c20);}
 { __m256 a3=_mm256_broadcast_ss(a+3);c30=_mm256_fmadd_ps(a3,b0,c30);}
 { __m256 a4=_mm256_broadcast_ss(a+4);c40=_mm256_fmadd_ps(a4,b0,c40);}
 { __m256 a5=_mm256_broadcast_ss(a+5);c50=_mm256_fmadd_ps(a5,b0,c50);}
 { __m256 a6=_mm256_broadcast_ss(a+6);c60=_mm256_fmadd_ps(a6,b0,c60);}
 { __m256 a7=_mm256_broadcast_ss(a+7);c70=_mm256_fmadd_ps(a7,b0,c70);}
 a+=MR;b+=NR;
 }
 _mm256_storeu_ps(c+0*stride+0,c00);
 _mm256_storeu_ps(c+1*stride+0,c10);
 _mm256_storeu_ps(c+2*stride+0,c20);
 _mm256_storeu_ps(c+3*stride+0,c30);
 _mm256_storeu_ps(c+4*stride+0,c40);
 _mm256_storeu_ps(c+5*stride+0,c50);
 _mm256_storeu_ps(c+6*stride+0,c60);
 _mm256_storeu_ps(c+7*stride+0,c70);
}
static void run(const float*A,const float*B,float*C){
 memset(C,0,(size_t)ROWS*N*4);
 float*ap=nullptr; if(posix_memalign((void**)&ap,4096,(size_t)ROWS*N*4))abort();
 constexpr int groups=ROWS/MR;
 for(int pc=0;pc<N;pc+=KC)for(int g=0;g<groups;++g)for(int k=0;k<KC;++k)for(int i=0;i<MR;++i)
  ap[(((size_t)(pc/KC)*groups+g)*KC+k)*MR+i]=A[(size_t)(g*MR+i)*N+pc+k];
 uint64_t begin=tick(),packing=0;
 for(int jc=0;jc<COLS;jc+=NC)for(int pc=0;pc<N;pc+=KC){
  uint64_t ps=tick();
  for(int q=0;q<NC;q+=NR)for(int k=0;k<KC;++k)for(int j=0;j<NR;j+=8)
   _mm256_store_ps(bp+((q/NR*KC+k)*NR)+j,_mm256_loadu_ps(B+(size_t)(pc+k)*N+jc+q+j));
  packing+=tick()-ps;
  for(int ic=0;ic<ROWS;ic+=MC)for(int q=0;q<NC;q+=NR)for(int r=0;r<MC;r+=MR)
   micro(ap+((size_t)(pc/KC)*groups+(ic+r)/MR)*KC*MR,bp+q*KC,C+(size_t)(ic+r)*N+jc+q,N);
 }
 uint64_t total=tick()-begin;
 for(int t=0;t<4;++t){int i=t*91,j=t*1999;double expected=0;
  for(int k=0;k<N;++k)expected+=(double)A[(size_t)i*N+k]*B[(size_t)k*N+j];
  if(fabs(expected-C[(size_t)i*N+j])>0.1)abort();
 }
 printf("m8n8u1 total_ticks=%llu pack_ticks=%llu core_ticks=%llu\n",(unsigned long long)total,(unsigned long long)packing,(unsigned long long)(total-packing));
 free(ap);
}
}
namespace m8n8u4 {
constexpr int MR=8,NR=8;
static inline void micro(const float*a,const float*b,float*c,int stride){
 __m256 c00=_mm256_loadu_ps(c+0*stride+0);
 __m256 c10=_mm256_loadu_ps(c+1*stride+0);
 __m256 c20=_mm256_loadu_ps(c+2*stride+0);
 __m256 c30=_mm256_loadu_ps(c+3*stride+0);
 __m256 c40=_mm256_loadu_ps(c+4*stride+0);
 __m256 c50=_mm256_loadu_ps(c+5*stride+0);
 __m256 c60=_mm256_loadu_ps(c+6*stride+0);
 __m256 c70=_mm256_loadu_ps(c+7*stride+0);
#pragma GCC unroll 4
 for(int k=0;k<KC;++k){
 __m256 b0=_mm256_load_ps(b+0);
 { __m256 a0=_mm256_broadcast_ss(a+0);c00=_mm256_fmadd_ps(a0,b0,c00);}
 { __m256 a1=_mm256_broadcast_ss(a+1);c10=_mm256_fmadd_ps(a1,b0,c10);}
 { __m256 a2=_mm256_broadcast_ss(a+2);c20=_mm256_fmadd_ps(a2,b0,c20);}
 { __m256 a3=_mm256_broadcast_ss(a+3);c30=_mm256_fmadd_ps(a3,b0,c30);}
 { __m256 a4=_mm256_broadcast_ss(a+4);c40=_mm256_fmadd_ps(a4,b0,c40);}
 { __m256 a5=_mm256_broadcast_ss(a+5);c50=_mm256_fmadd_ps(a5,b0,c50);}
 { __m256 a6=_mm256_broadcast_ss(a+6);c60=_mm256_fmadd_ps(a6,b0,c60);}
 { __m256 a7=_mm256_broadcast_ss(a+7);c70=_mm256_fmadd_ps(a7,b0,c70);}
 a+=MR;b+=NR;
 }
 _mm256_storeu_ps(c+0*stride+0,c00);
 _mm256_storeu_ps(c+1*stride+0,c10);
 _mm256_storeu_ps(c+2*stride+0,c20);
 _mm256_storeu_ps(c+3*stride+0,c30);
 _mm256_storeu_ps(c+4*stride+0,c40);
 _mm256_storeu_ps(c+5*stride+0,c50);
 _mm256_storeu_ps(c+6*stride+0,c60);
 _mm256_storeu_ps(c+7*stride+0,c70);
}
static void run(const float*A,const float*B,float*C){
 memset(C,0,(size_t)ROWS*N*4);
 float*ap=nullptr; if(posix_memalign((void**)&ap,4096,(size_t)ROWS*N*4))abort();
 constexpr int groups=ROWS/MR;
 for(int pc=0;pc<N;pc+=KC)for(int g=0;g<groups;++g)for(int k=0;k<KC;++k)for(int i=0;i<MR;++i)
  ap[(((size_t)(pc/KC)*groups+g)*KC+k)*MR+i]=A[(size_t)(g*MR+i)*N+pc+k];
 uint64_t begin=tick(),packing=0;
 for(int jc=0;jc<COLS;jc+=NC)for(int pc=0;pc<N;pc+=KC){
  uint64_t ps=tick();
  for(int q=0;q<NC;q+=NR)for(int k=0;k<KC;++k)for(int j=0;j<NR;j+=8)
   _mm256_store_ps(bp+((q/NR*KC+k)*NR)+j,_mm256_loadu_ps(B+(size_t)(pc+k)*N+jc+q+j));
  packing+=tick()-ps;
  for(int ic=0;ic<ROWS;ic+=MC)for(int q=0;q<NC;q+=NR)for(int r=0;r<MC;r+=MR)
   micro(ap+((size_t)(pc/KC)*groups+(ic+r)/MR)*KC*MR,bp+q*KC,C+(size_t)(ic+r)*N+jc+q,N);
 }
 uint64_t total=tick()-begin;
 for(int t=0;t<4;++t){int i=t*91,j=t*1999;double expected=0;
  for(int k=0;k<N;++k)expected+=(double)A[(size_t)i*N+k]*B[(size_t)k*N+j];
  if(fabs(expected-C[(size_t)i*N+j])>0.1)abort();
 }
 printf("m8n8u4 total_ticks=%llu pack_ticks=%llu core_ticks=%llu\n",(unsigned long long)total,(unsigned long long)packing,(unsigned long long)(total-packing));
 free(ap);
}
}
namespace m4n16u1 {
constexpr int MR=4,NR=16;
static inline void micro(const float*a,const float*b,float*c,int stride){
 __m256 c00=_mm256_loadu_ps(c+0*stride+0);
 __m256 c01=_mm256_loadu_ps(c+0*stride+8);
 __m256 c10=_mm256_loadu_ps(c+1*stride+0);
 __m256 c11=_mm256_loadu_ps(c+1*stride+8);
 __m256 c20=_mm256_loadu_ps(c+2*stride+0);
 __m256 c21=_mm256_loadu_ps(c+2*stride+8);
 __m256 c30=_mm256_loadu_ps(c+3*stride+0);
 __m256 c31=_mm256_loadu_ps(c+3*stride+8);
#pragma GCC unroll 1
 for(int k=0;k<KC;++k){
 __m256 b0=_mm256_load_ps(b+0);
 __m256 b1=_mm256_load_ps(b+8);
 { __m256 a0=_mm256_broadcast_ss(a+0);c00=_mm256_fmadd_ps(a0,b0,c00);c01=_mm256_fmadd_ps(a0,b1,c01);}
 { __m256 a1=_mm256_broadcast_ss(a+1);c10=_mm256_fmadd_ps(a1,b0,c10);c11=_mm256_fmadd_ps(a1,b1,c11);}
 { __m256 a2=_mm256_broadcast_ss(a+2);c20=_mm256_fmadd_ps(a2,b0,c20);c21=_mm256_fmadd_ps(a2,b1,c21);}
 { __m256 a3=_mm256_broadcast_ss(a+3);c30=_mm256_fmadd_ps(a3,b0,c30);c31=_mm256_fmadd_ps(a3,b1,c31);}
 a+=MR;b+=NR;
 }
 _mm256_storeu_ps(c+0*stride+0,c00);
 _mm256_storeu_ps(c+0*stride+8,c01);
 _mm256_storeu_ps(c+1*stride+0,c10);
 _mm256_storeu_ps(c+1*stride+8,c11);
 _mm256_storeu_ps(c+2*stride+0,c20);
 _mm256_storeu_ps(c+2*stride+8,c21);
 _mm256_storeu_ps(c+3*stride+0,c30);
 _mm256_storeu_ps(c+3*stride+8,c31);
}
static void run(const float*A,const float*B,float*C){
 memset(C,0,(size_t)ROWS*N*4);
 float*ap=nullptr; if(posix_memalign((void**)&ap,4096,(size_t)ROWS*N*4))abort();
 constexpr int groups=ROWS/MR;
 for(int pc=0;pc<N;pc+=KC)for(int g=0;g<groups;++g)for(int k=0;k<KC;++k)for(int i=0;i<MR;++i)
  ap[(((size_t)(pc/KC)*groups+g)*KC+k)*MR+i]=A[(size_t)(g*MR+i)*N+pc+k];
 uint64_t begin=tick(),packing=0;
 for(int jc=0;jc<COLS;jc+=NC)for(int pc=0;pc<N;pc+=KC){
  uint64_t ps=tick();
  for(int q=0;q<NC;q+=NR)for(int k=0;k<KC;++k)for(int j=0;j<NR;j+=8)
   _mm256_store_ps(bp+((q/NR*KC+k)*NR)+j,_mm256_loadu_ps(B+(size_t)(pc+k)*N+jc+q+j));
  packing+=tick()-ps;
  for(int ic=0;ic<ROWS;ic+=MC)for(int q=0;q<NC;q+=NR)for(int r=0;r<MC;r+=MR)
   micro(ap+((size_t)(pc/KC)*groups+(ic+r)/MR)*KC*MR,bp+q*KC,C+(size_t)(ic+r)*N+jc+q,N);
 }
 uint64_t total=tick()-begin;
 for(int t=0;t<4;++t){int i=t*91,j=t*1999;double expected=0;
  for(int k=0;k<N;++k)expected+=(double)A[(size_t)i*N+k]*B[(size_t)k*N+j];
  if(fabs(expected-C[(size_t)i*N+j])>0.1)abort();
 }
 printf("m4n16u1 total_ticks=%llu pack_ticks=%llu core_ticks=%llu\n",(unsigned long long)total,(unsigned long long)packing,(unsigned long long)(total-packing));
 free(ap);
}
}
namespace m4n16u2 {
constexpr int MR=4,NR=16;
static inline void micro(const float*a,const float*b,float*c,int stride){
 __m256 c00=_mm256_loadu_ps(c+0*stride+0);
 __m256 c01=_mm256_loadu_ps(c+0*stride+8);
 __m256 c10=_mm256_loadu_ps(c+1*stride+0);
 __m256 c11=_mm256_loadu_ps(c+1*stride+8);
 __m256 c20=_mm256_loadu_ps(c+2*stride+0);
 __m256 c21=_mm256_loadu_ps(c+2*stride+8);
 __m256 c30=_mm256_loadu_ps(c+3*stride+0);
 __m256 c31=_mm256_loadu_ps(c+3*stride+8);
#pragma GCC unroll 2
 for(int k=0;k<KC;++k){
 __m256 b0=_mm256_load_ps(b+0);
 __m256 b1=_mm256_load_ps(b+8);
 { __m256 a0=_mm256_broadcast_ss(a+0);c00=_mm256_fmadd_ps(a0,b0,c00);c01=_mm256_fmadd_ps(a0,b1,c01);}
 { __m256 a1=_mm256_broadcast_ss(a+1);c10=_mm256_fmadd_ps(a1,b0,c10);c11=_mm256_fmadd_ps(a1,b1,c11);}
 { __m256 a2=_mm256_broadcast_ss(a+2);c20=_mm256_fmadd_ps(a2,b0,c20);c21=_mm256_fmadd_ps(a2,b1,c21);}
 { __m256 a3=_mm256_broadcast_ss(a+3);c30=_mm256_fmadd_ps(a3,b0,c30);c31=_mm256_fmadd_ps(a3,b1,c31);}
 a+=MR;b+=NR;
 }
 _mm256_storeu_ps(c+0*stride+0,c00);
 _mm256_storeu_ps(c+0*stride+8,c01);
 _mm256_storeu_ps(c+1*stride+0,c10);
 _mm256_storeu_ps(c+1*stride+8,c11);
 _mm256_storeu_ps(c+2*stride+0,c20);
 _mm256_storeu_ps(c+2*stride+8,c21);
 _mm256_storeu_ps(c+3*stride+0,c30);
 _mm256_storeu_ps(c+3*stride+8,c31);
}
static void run(const float*A,const float*B,float*C){
 memset(C,0,(size_t)ROWS*N*4);
 float*ap=nullptr; if(posix_memalign((void**)&ap,4096,(size_t)ROWS*N*4))abort();
 constexpr int groups=ROWS/MR;
 for(int pc=0;pc<N;pc+=KC)for(int g=0;g<groups;++g)for(int k=0;k<KC;++k)for(int i=0;i<MR;++i)
  ap[(((size_t)(pc/KC)*groups+g)*KC+k)*MR+i]=A[(size_t)(g*MR+i)*N+pc+k];
 uint64_t begin=tick(),packing=0;
 for(int jc=0;jc<COLS;jc+=NC)for(int pc=0;pc<N;pc+=KC){
  uint64_t ps=tick();
  for(int q=0;q<NC;q+=NR)for(int k=0;k<KC;++k)for(int j=0;j<NR;j+=8)
   _mm256_store_ps(bp+((q/NR*KC+k)*NR)+j,_mm256_loadu_ps(B+(size_t)(pc+k)*N+jc+q+j));
  packing+=tick()-ps;
  for(int ic=0;ic<ROWS;ic+=MC)for(int q=0;q<NC;q+=NR)for(int r=0;r<MC;r+=MR)
   micro(ap+((size_t)(pc/KC)*groups+(ic+r)/MR)*KC*MR,bp+q*KC,C+(size_t)(ic+r)*N+jc+q,N);
 }
 uint64_t total=tick()-begin;
 for(int t=0;t<4;++t){int i=t*91,j=t*1999;double expected=0;
  for(int k=0;k<N;++k)expected+=(double)A[(size_t)i*N+k]*B[(size_t)k*N+j];
  if(fabs(expected-C[(size_t)i*N+j])>0.1)abort();
 }
 printf("m4n16u2 total_ticks=%llu pack_ticks=%llu core_ticks=%llu\n",(unsigned long long)total,(unsigned long long)packing,(unsigned long long)(total-packing));
 free(ap);
}
}
namespace m4n16u4 {
constexpr int MR=4,NR=16;
static inline void micro(const float*a,const float*b,float*c,int stride){
 __m256 c00=_mm256_loadu_ps(c+0*stride+0);
 __m256 c01=_mm256_loadu_ps(c+0*stride+8);
 __m256 c10=_mm256_loadu_ps(c+1*stride+0);
 __m256 c11=_mm256_loadu_ps(c+1*stride+8);
 __m256 c20=_mm256_loadu_ps(c+2*stride+0);
 __m256 c21=_mm256_loadu_ps(c+2*stride+8);
 __m256 c30=_mm256_loadu_ps(c+3*stride+0);
 __m256 c31=_mm256_loadu_ps(c+3*stride+8);
#pragma GCC unroll 4
 for(int k=0;k<KC;++k){
 __m256 b0=_mm256_load_ps(b+0);
 __m256 b1=_mm256_load_ps(b+8);
 { __m256 a0=_mm256_broadcast_ss(a+0);c00=_mm256_fmadd_ps(a0,b0,c00);c01=_mm256_fmadd_ps(a0,b1,c01);}
 { __m256 a1=_mm256_broadcast_ss(a+1);c10=_mm256_fmadd_ps(a1,b0,c10);c11=_mm256_fmadd_ps(a1,b1,c11);}
 { __m256 a2=_mm256_broadcast_ss(a+2);c20=_mm256_fmadd_ps(a2,b0,c20);c21=_mm256_fmadd_ps(a2,b1,c21);}
 { __m256 a3=_mm256_broadcast_ss(a+3);c30=_mm256_fmadd_ps(a3,b0,c30);c31=_mm256_fmadd_ps(a3,b1,c31);}
 a+=MR;b+=NR;
 }
 _mm256_storeu_ps(c+0*stride+0,c00);
 _mm256_storeu_ps(c+0*stride+8,c01);
 _mm256_storeu_ps(c+1*stride+0,c10);
 _mm256_storeu_ps(c+1*stride+8,c11);
 _mm256_storeu_ps(c+2*stride+0,c20);
 _mm256_storeu_ps(c+2*stride+8,c21);
 _mm256_storeu_ps(c+3*stride+0,c30);
 _mm256_storeu_ps(c+3*stride+8,c31);
}
static void run(const float*A,const float*B,float*C){
 memset(C,0,(size_t)ROWS*N*4);
 float*ap=nullptr; if(posix_memalign((void**)&ap,4096,(size_t)ROWS*N*4))abort();
 constexpr int groups=ROWS/MR;
 for(int pc=0;pc<N;pc+=KC)for(int g=0;g<groups;++g)for(int k=0;k<KC;++k)for(int i=0;i<MR;++i)
  ap[(((size_t)(pc/KC)*groups+g)*KC+k)*MR+i]=A[(size_t)(g*MR+i)*N+pc+k];
 uint64_t begin=tick(),packing=0;
 for(int jc=0;jc<COLS;jc+=NC)for(int pc=0;pc<N;pc+=KC){
  uint64_t ps=tick();
  for(int q=0;q<NC;q+=NR)for(int k=0;k<KC;++k)for(int j=0;j<NR;j+=8)
   _mm256_store_ps(bp+((q/NR*KC+k)*NR)+j,_mm256_loadu_ps(B+(size_t)(pc+k)*N+jc+q+j));
  packing+=tick()-ps;
  for(int ic=0;ic<ROWS;ic+=MC)for(int q=0;q<NC;q+=NR)for(int r=0;r<MC;r+=MR)
   micro(ap+((size_t)(pc/KC)*groups+(ic+r)/MR)*KC*MR,bp+q*KC,C+(size_t)(ic+r)*N+jc+q,N);
 }
 uint64_t total=tick()-begin;
 for(int t=0;t<4;++t){int i=t*91,j=t*1999;double expected=0;
  for(int k=0;k<N;++k)expected+=(double)A[(size_t)i*N+k]*B[(size_t)k*N+j];
  if(fabs(expected-C[(size_t)i*N+j])>0.1)abort();
 }
 printf("m4n16u4 total_ticks=%llu pack_ticks=%llu core_ticks=%llu\n",(unsigned long long)total,(unsigned long long)packing,(unsigned long long)(total-packing));
 free(ap);
}
}
namespace m4n16u8 {
constexpr int MR=4,NR=16;
static inline void micro(const float*a,const float*b,float*c,int stride){
 __m256 c00=_mm256_loadu_ps(c+0*stride+0);
 __m256 c01=_mm256_loadu_ps(c+0*stride+8);
 __m256 c10=_mm256_loadu_ps(c+1*stride+0);
 __m256 c11=_mm256_loadu_ps(c+1*stride+8);
 __m256 c20=_mm256_loadu_ps(c+2*stride+0);
 __m256 c21=_mm256_loadu_ps(c+2*stride+8);
 __m256 c30=_mm256_loadu_ps(c+3*stride+0);
 __m256 c31=_mm256_loadu_ps(c+3*stride+8);
#pragma GCC unroll 8
 for(int k=0;k<KC;++k){
 __m256 b0=_mm256_load_ps(b+0);
 __m256 b1=_mm256_load_ps(b+8);
 { __m256 a0=_mm256_broadcast_ss(a+0);c00=_mm256_fmadd_ps(a0,b0,c00);c01=_mm256_fmadd_ps(a0,b1,c01);}
 { __m256 a1=_mm256_broadcast_ss(a+1);c10=_mm256_fmadd_ps(a1,b0,c10);c11=_mm256_fmadd_ps(a1,b1,c11);}
 { __m256 a2=_mm256_broadcast_ss(a+2);c20=_mm256_fmadd_ps(a2,b0,c20);c21=_mm256_fmadd_ps(a2,b1,c21);}
 { __m256 a3=_mm256_broadcast_ss(a+3);c30=_mm256_fmadd_ps(a3,b0,c30);c31=_mm256_fmadd_ps(a3,b1,c31);}
 a+=MR;b+=NR;
 }
 _mm256_storeu_ps(c+0*stride+0,c00);
 _mm256_storeu_ps(c+0*stride+8,c01);
 _mm256_storeu_ps(c+1*stride+0,c10);
 _mm256_storeu_ps(c+1*stride+8,c11);
 _mm256_storeu_ps(c+2*stride+0,c20);
 _mm256_storeu_ps(c+2*stride+8,c21);
 _mm256_storeu_ps(c+3*stride+0,c30);
 _mm256_storeu_ps(c+3*stride+8,c31);
}
static void run(const float*A,const float*B,float*C){
 memset(C,0,(size_t)ROWS*N*4);
 float*ap=nullptr; if(posix_memalign((void**)&ap,4096,(size_t)ROWS*N*4))abort();
 constexpr int groups=ROWS/MR;
 for(int pc=0;pc<N;pc+=KC)for(int g=0;g<groups;++g)for(int k=0;k<KC;++k)for(int i=0;i<MR;++i)
  ap[(((size_t)(pc/KC)*groups+g)*KC+k)*MR+i]=A[(size_t)(g*MR+i)*N+pc+k];
 uint64_t begin=tick(),packing=0;
 for(int jc=0;jc<COLS;jc+=NC)for(int pc=0;pc<N;pc+=KC){
  uint64_t ps=tick();
  for(int q=0;q<NC;q+=NR)for(int k=0;k<KC;++k)for(int j=0;j<NR;j+=8)
   _mm256_store_ps(bp+((q/NR*KC+k)*NR)+j,_mm256_loadu_ps(B+(size_t)(pc+k)*N+jc+q+j));
  packing+=tick()-ps;
  for(int ic=0;ic<ROWS;ic+=MC)for(int q=0;q<NC;q+=NR)for(int r=0;r<MC;r+=MR)
   micro(ap+((size_t)(pc/KC)*groups+(ic+r)/MR)*KC*MR,bp+q*KC,C+(size_t)(ic+r)*N+jc+q,N);
 }
 uint64_t total=tick()-begin;
 for(int t=0;t<4;++t){int i=t*91,j=t*1999;double expected=0;
  for(int k=0;k<N;++k)expected+=(double)A[(size_t)i*N+k]*B[(size_t)k*N+j];
  if(fabs(expected-C[(size_t)i*N+j])>0.1)abort();
 }
 printf("m4n16u8 total_ticks=%llu pack_ticks=%llu core_ticks=%llu\n",(unsigned long long)total,(unsigned long long)packing,(unsigned long long)(total-packing));
 free(ap);
}
}
int main(){
 float *A=nullptr,*B=nullptr,*C=nullptr;
 if(posix_memalign((void**)&A,4096,(size_t)ROWS*N*4)||posix_memalign((void**)&B,4096,(size_t)N*N*4)||posix_memalign((void**)&C,4096,(size_t)ROWS*N*4))abort();
 for(size_t i=0;i<(size_t)ROWS*N;++i)A[i]=(float)((int)((i*17+43)%257)-128)*(1.0f/128);
 for(size_t i=0;i<(size_t)N*N;++i)B[i]=(float)((int)((i*23+71)%263)-131)*(1.0f/131);
m6n16u1::run(A,B,C);
m6n16u2::run(A,B,C);
m6n16u4::run(A,B,C);
m6n16u8::run(A,B,C);
m4n24u1::run(A,B,C);
m4n24u2::run(A,B,C);
m4n24u4::run(A,B,C);
m8n8u1::run(A,B,C);
m8n8u4::run(A,B,C);
m4n16u1::run(A,B,C);
m4n16u2::run(A,B,C);
m4n16u4::run(A,B,C);
m4n16u8::run(A,B,C);
free(C);free(B);free(A);return 0;
}

评测结果