代码

/*
Credits:
- saffah_cc_v41_260924, https://duck.ac/submission/99846:
  fused leaf operand sums/packing and one-pass four-quadrant Strassen combine.
- saffah_cc_v41_agg1, https://duck.ac/submission/99875 and
  https://duck.ac/submission/99950: non-temporal intermediate/output stores.
- saffah_cc_v41_260924, https://duck.ac/submission/99141, and
  saffah_cc_v41_agg1, https://duck.ac/submission/99167:
  double-precision 6x8 blocking, explicit 256-bit memory instructions, and
  row-contiguous A packing and row-first B packing were consulted for this problem's baseline.
- pdoom, https://duck.ac/submission/99966: preceding packing/assembly framework.
The double microkernel, output-backed temporary placement, and hybrid
Strassen/Winograd recursion below were implemented here.
No quantization is used: input, arithmetic, and output remain binary64.
*/
#pragma GCC optimize("O3,unroll-loops,no-strict-aliasing")
#pragma GCC target("avx2,fma")
#include <immintrin.h>
#include <stddef.h>
#ifndef BASE
#define BASE 512
#endif
#ifndef MC
#define MC 48
#endif
#ifndef NTMIN
#define NTMIN 512
#endif
#ifndef SUMNT
#define SUMNT 1024
#endif
static constexpr int MR=6,NR=8;
alignas(4096) static double workspace[27*1024*1024];
alignas(64) static double ap[(MC+MR)*BASE],bp[BASE*BASE];
struct Mat { const double *x; int ld; const double *y=nullptr; int ly=0; int sign=0; };
static inline __m256d ld(Mat a,int i,int j) {
    __m256d x=_mm256_load_pd(a.x+(size_t)i*a.ld+j);
    if(a.sign) {
        __m256d y=_mm256_load_pd(a.y+(size_t)i*a.ly+j);
        x=a.sign>0?_mm256_add_pd(x,y):_mm256_sub_pd(x,y);
    }
    return x;
}
static inline double scalar(Mat a,int i,int k) {
    double x=a.x[(size_t)i*a.ld+k];
    if(a.sign)x+=a.sign>0?a.y[(size_t)i*a.ly+k]:-a.y[(size_t)i*a.ly+k];
    return x;
}
static __attribute__((noinline)) void kernel(int n,const double *a,const double *b,double *c,int ldc,int rows) {
    alignas(32) double tmp[48];
    double *out=rows==MR?c:tmp;
    ptrdiff_t stride=rows==MR?ldc*8:64;
    int k=n;
    asm volatile (
        "vxorpd %%ymm0, %%ymm0, %%ymm0\n\t"
        "vxorpd %%ymm1, %%ymm1, %%ymm1\n\t"
        "vxorpd %%ymm2, %%ymm2, %%ymm2\n\t"
        "vxorpd %%ymm3, %%ymm3, %%ymm3\n\t"
        "vxorpd %%ymm4, %%ymm4, %%ymm4\n\t"
        "vxorpd %%ymm5, %%ymm5, %%ymm5\n\t"
        "vxorpd %%ymm6, %%ymm6, %%ymm6\n\t"
        "vxorpd %%ymm7, %%ymm7, %%ymm7\n\t"
        "vxorpd %%ymm8, %%ymm8, %%ymm8\n\t"
        "vxorpd %%ymm9, %%ymm9, %%ymm9\n\t"
        "vxorpd %%ymm10, %%ymm10, %%ymm10\n\t"
        "vxorpd %%ymm11, %%ymm11, %%ymm11\n\t"
        ".p2align 5\n\t"
        "1:\n\t"
        "vmovapd 0(%[b]), %%ymm12\n\t"
        "vmovapd 32(%[b]), %%ymm13\n\t"
        "vbroadcastsd %c[o0]+0(%[a]), %%ymm14\n\t"
        "vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
        "vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
        "vbroadcastsd %c[o1]+0(%[a]), %%ymm14\n\t"
        "vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
        "vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
        "vbroadcastsd %c[o2]+0(%[a]), %%ymm14\n\t"
        "vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
        "vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
        "vbroadcastsd %c[o3]+0(%[a]), %%ymm14\n\t"
        "vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
        "vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
        "vbroadcastsd %c[o4]+0(%[a]), %%ymm14\n\t"
        "vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
        "vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
        "vbroadcastsd %c[o5]+0(%[a]), %%ymm14\n\t"
        "vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
        "vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
        "vmovapd 64(%[b]), %%ymm12\n\t"
        "vmovapd 96(%[b]), %%ymm13\n\t"
        "vbroadcastsd %c[o0]+8(%[a]), %%ymm14\n\t"
        "vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
        "vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
        "vbroadcastsd %c[o1]+8(%[a]), %%ymm14\n\t"
        "vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
        "vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
        "vbroadcastsd %c[o2]+8(%[a]), %%ymm14\n\t"
        "vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
        "vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
        "vbroadcastsd %c[o3]+8(%[a]), %%ymm14\n\t"
        "vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
        "vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
        "vbroadcastsd %c[o4]+8(%[a]), %%ymm14\n\t"
        "vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
        "vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
        "vbroadcastsd %c[o5]+8(%[a]), %%ymm14\n\t"
        "vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
        "vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
        "vmovapd 128(%[b]), %%ymm12\n\t"
        "vmovapd 160(%[b]), %%ymm13\n\t"
        "vbroadcastsd %c[o0]+16(%[a]), %%ymm14\n\t"
        "vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
        "vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
        "vbroadcastsd %c[o1]+16(%[a]), %%ymm14\n\t"
        "vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
        "vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
        "vbroadcastsd %c[o2]+16(%[a]), %%ymm14\n\t"
        "vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
        "vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
        "vbroadcastsd %c[o3]+16(%[a]), %%ymm14\n\t"
        "vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
        "vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
        "vbroadcastsd %c[o4]+16(%[a]), %%ymm14\n\t"
        "vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
        "vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
        "vbroadcastsd %c[o5]+16(%[a]), %%ymm14\n\t"
        "vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
        "vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
        "vmovapd 192(%[b]), %%ymm12\n\t"
        "vmovapd 224(%[b]), %%ymm13\n\t"
        "vbroadcastsd %c[o0]+24(%[a]), %%ymm14\n\t"
        "vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
        "vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
        "vbroadcastsd %c[o1]+24(%[a]), %%ymm14\n\t"
        "vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
        "vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
        "vbroadcastsd %c[o2]+24(%[a]), %%ymm14\n\t"
        "vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
        "vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
        "vbroadcastsd %c[o3]+24(%[a]), %%ymm14\n\t"
        "vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
        "vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
        "vbroadcastsd %c[o4]+24(%[a]), %%ymm14\n\t"
        "vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
        "vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
        "vbroadcastsd %c[o5]+24(%[a]), %%ymm14\n\t"
        "vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
        "vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
        "add $32, %[a]\n\t"
        "add $256, %[b]\n\t"
        "sub $4, %[k]\n\t"
        "jnz 1b\n\t"
        "test %[nt], %[nt]\n\t"
        "jz 2f\n\t"
        "vmovntpd %%ymm0, (%[out])\n\t"
        "vmovntpd %%ymm1, 32(%[out])\n\t"
        "add %[stride], %[out]\n\t"
        "vmovntpd %%ymm2, (%[out])\n\t"
        "vmovntpd %%ymm3, 32(%[out])\n\t"
        "add %[stride], %[out]\n\t"
        "vmovntpd %%ymm4, (%[out])\n\t"
        "vmovntpd %%ymm5, 32(%[out])\n\t"
        "add %[stride], %[out]\n\t"
        "vmovntpd %%ymm6, (%[out])\n\t"
        "vmovntpd %%ymm7, 32(%[out])\n\t"
        "add %[stride], %[out]\n\t"
        "vmovntpd %%ymm8, (%[out])\n\t"
        "vmovntpd %%ymm9, 32(%[out])\n\t"
        "add %[stride], %[out]\n\t"
        "vmovntpd %%ymm10, (%[out])\n\t"
        "vmovntpd %%ymm11, 32(%[out])\n\t"
        "jmp 3f\n\t"
        "2:\n\t"
        "vmovapd %%ymm0, (%[out])\n\t"
        "vmovapd %%ymm1, 32(%[out])\n\t"
        "add %[stride], %[out]\n\t"
        "vmovapd %%ymm2, (%[out])\n\t"
        "vmovapd %%ymm3, 32(%[out])\n\t"
        "add %[stride], %[out]\n\t"
        "vmovapd %%ymm4, (%[out])\n\t"
        "vmovapd %%ymm5, 32(%[out])\n\t"
        "add %[stride], %[out]\n\t"
        "vmovapd %%ymm6, (%[out])\n\t"
        "vmovapd %%ymm7, 32(%[out])\n\t"
        "add %[stride], %[out]\n\t"
        "vmovapd %%ymm8, (%[out])\n\t"
        "vmovapd %%ymm9, 32(%[out])\n\t"
        "add %[stride], %[out]\n\t"
        "vmovapd %%ymm10, (%[out])\n\t"
        "vmovapd %%ymm11, 32(%[out])\n\t"
        "3:\n\t"
        : [a] "+&r"(a), [b] "+&r"(b), [out] "+&r"(out), [k] "+&r"(k)
        : [stride] "r"(stride), [nt] "r"(int(n>=NTMIN&&rows==MR)), [o0] "i"(0*BASE*8), [o1] "i"(1*BASE*8), [o2] "i"(2*BASE*8), [o3] "i"(3*BASE*8), [o4] "i"(4*BASE*8), [o5] "i"(5*BASE*8)
        : "memory", "cc", "ymm0","ymm1","ymm2","ymm3","ymm4","ymm5","ymm6","ymm7","ymm8","ymm9","ymm10","ymm11","ymm12","ymm13","ymm14","ymm15");
    if(rows!=MR)for(int i=0;i<rows;i++)for(int j=0;j<8;j+=4)
        _mm256_store_pd(c+(size_t)i*ldc+j,_mm256_load_pd(tmp+i*8+j));
}
static void leaf(int n,Mat a,Mat b,double *c,int ldc) {
    for(int k=0;k<n;k++)for(int j=0;j<n;j+=8) {
        double *p=bp+(size_t)j*n+k*8;
        _mm256_store_pd(p,ld(b,k,j));_mm256_store_pd(p+4,ld(b,k,j+4));
    }
    for(int ic=0;ic<n;ic+=MC) {
        int mc=n-ic<MC?n-ic:MC;
        for(int r=0;r<mc;r+=MR) {
            int rows=mc-r<MR?mc-r:MR;
            for(int i=0;i<rows;i++)for(int k=0;k<n;k+=4)
                _mm256_store_pd(ap+(size_t)(r+i)*n+k,ld(a,ic+r+i,k));
            for(int i=rows;i<MR;i++)for(int k=0;k<n;k+=4)
                _mm256_store_pd(ap+(size_t)(r+i)*n+k,_mm256_setzero_pd());
        }
        for(int j=0;j<n;j+=NR)for(int r=0;r<mc;r+=MR)
            kernel(n,ap+(size_t)r*n,bp+(size_t)j*n,c+(size_t)(ic+r)*ldc+j,ldc,mc-r<MR?mc-r:MR);
    }
    if(n>=NTMIN)_mm_sfence();
}
static void sum(int h,double *dst,Mat a,Mat b,int sign) {
    for(int i=0;i<h;i++)for(int j=0;j<h;j+=4) {
        __m256d v=sign>0?_mm256_add_pd(ld(a,i,j),ld(b,i,j)):_mm256_sub_pd(ld(a,i,j),ld(b,i,j));
        if(h>=SUMNT)_mm256_stream_pd(dst+(size_t)i*h+j,v);
        else _mm256_store_pd(dst+(size_t)i*h+j,v);
    }
    if(h>=SUMNT)_mm_sfence();
}
static void combine(int h,double *c,int ldc,const double *m,bool wino=false) {
    size_t z=(size_t)h*h;
    for(int i=0;i<h;i++)for(int j=0;j<h;j+=4) {
        size_t t=(size_t)i*h+j;
        double *c11=c+(size_t)i*ldc+j,*c12=c11+h,*c21=c+(size_t)(i+h)*ldc+j,*c22=c21+h;
        __m256d m1=_mm256_load_pd(c11),m2=_mm256_load_pd(c21),m3=_mm256_load_pd(c12),m4=_mm256_load_pd(c22),
          m5=_mm256_load_pd(m+t),m6=_mm256_load_pd(m+z+t),m7=_mm256_load_pd(m+2*z+t);
        __m256d r11,r12,r21,r22;
        if(wino) {
            __m256d u2=_mm256_add_pd(m1,m6),u3=_mm256_add_pd(u2,m7),u4=_mm256_add_pd(u2,m5);
            r11=_mm256_add_pd(m1,m2);r12=_mm256_add_pd(u4,m3);r21=_mm256_sub_pd(u3,m4);r22=_mm256_add_pd(u3,m5);
        } else {
            r11=_mm256_add_pd(_mm256_sub_pd(_mm256_add_pd(m1,m4),m5),m7);r12=_mm256_add_pd(m3,m5);r21=_mm256_add_pd(m2,m4);
            r22=_mm256_add_pd(_mm256_add_pd(_mm256_sub_pd(m1,m2),m3),m6);
        }
        _mm256_store_pd(c11,r11);_mm256_store_pd(c12,r12);_mm256_store_pd(c21,r21);_mm256_store_pd(c22,r22);
    }
}
static void rec(int n,Mat a,Mat b,double *c,int ldc,double *work) {
    if(n<=BASE){leaf(n,a,b,c,ldc);return;}
    int h=n/2;size_t z=(size_t)h*h;
    Mat aa[4]={{a.x,a.ld},{a.x+h,a.ld},{a.x+(size_t)h*a.ld,a.ld},{a.x+(size_t)h*a.ld+h,a.ld}};
    Mat bb[4]={{b.x,b.ld},{b.x+h,b.ld},{b.x+(size_t)h*b.ld,b.ld},{b.x+(size_t)h*b.ld+h,b.ld}};
    double *m=work,*p=m+3*z,*q=p+z,*sub=q+z;
    double *out[7]={c,c+(size_t)h*ldc,c+h,c+(size_t)h*ldc+h,m,m+z,m+2*z};
    if(n>BASE*2) {
        rec(h,aa[0],bb[0],out[0],ldc,sub);
        rec(h,aa[1],bb[2],out[1],ldc,sub);
        sum(h,p,aa[2],aa[3],1);sum(h,q,bb[1],bb[0],-1);
        rec(h,{p,h},{q,h},out[4],h,sub);
        sum(h,p,{p,h},aa[0],-1);sum(h,q,bb[3],{q,h},-1);
        rec(h,{p,h},{q,h},out[5],h,sub);
        sum(h,p,aa[1],{p,h},-1);
        rec(h,{p,h},bb[3],out[2],ldc,sub);
        sum(h,q,{q,h},bb[2],-1);
        rec(h,aa[3],{q,h},out[3],ldc,sub);
        sum(h,p,aa[0],aa[2],-1);sum(h,q,bb[3],bb[1],-1);
        rec(h,{p,h},{q,h},out[6],h,sub);
        combine(h,c,ldc,m,true);
        return;
    }
    const int ai[7]={0,2,0,3,0,2,1},aj[7]={3,3,-1,-1,1,0,3},as[7]={1,1,0,0,1,-1,-1};
    const int bi[7]={0,0,1,2,3,0,2},bj[7]={3,-1,3,0,-1,1,3},bs[7]={1,0,-1,-1,0,1,1};
    for(int r=0;r<7;r++) {
        Mat x=aa[ai[r]],y=bb[bi[r]];
        if(h<=BASE) {
            if(as[r]){x.y=aa[aj[r]].x;x.ly=a.ld;x.sign=as[r];}
            if(bs[r]){y.y=bb[bj[r]].x;y.ly=b.ld;y.sign=bs[r];}
            leaf(h,x,y,out[r],r<4?ldc:h);
        } else {
            if(as[r]){sum(h,p,x,aa[aj[r]],as[r]);x={p,h};}
            if(bs[r]){sum(h,q,y,bb[bj[r]],bs[r]);y={q,h};}
            rec(h,x,y,out[r],r<4?ldc:h,sub);
        }
    }
    combine(h,c,ldc,m);
}
void matrix_multiply(int n,const double *a,const double *b,double *c) {
    rec(n,{a,n},{b,n},c,n,workspace);
    _mm_sfence();
}

评测结果