代码

// DUCK_PUBLIC_SURROGATE_BENCHMARK
// Independent arithmetic matrix pattern; contains no production data or seeds.
#pragma GCC optimize("O3")
#pragma GCC target("avx2,fma")
#include <immintrin.h>
#include <stdint.h>
#include <stddef.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <stdio.h>
namespace check_float {
constexpr int N=8192,MR=4,NR=16,MC=128,NC=128,KC=256;
static float apack[MC*KC] __attribute__((aligned(4096)));
static float bpack[NC*KC] __attribute__((aligned(4096)));
static inline void micro(const float*a,const float*b,float*c,int stride) {
    __m256 c00=_mm256_loadu_ps(c),c01=_mm256_loadu_ps(c+8);
    __m256 c10=_mm256_loadu_ps(c+stride),c11=_mm256_loadu_ps(c+stride+8);
    __m256 c20=_mm256_loadu_ps(c+2*stride),c21=_mm256_loadu_ps(c+2*stride+8);
    __m256 c30=_mm256_loadu_ps(c+3*stride),c31=_mm256_loadu_ps(c+3*stride+8);
#pragma GCC unroll 8
    for(int k=0;k<KC;++k){
        __m256 b0=_mm256_loadu_ps(b),b1=_mm256_loadu_ps(b+8);
        __m256 a0=_mm256_broadcast_ss(a);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,c00);_mm256_storeu_ps(c+8,c01);
    _mm256_storeu_ps(c+stride,c10);_mm256_storeu_ps(c+stride+8,c11);
    _mm256_storeu_ps(c+2*stride,c20);_mm256_storeu_ps(c+2*stride+8,c21);
    _mm256_storeu_ps(c+3*stride,c30);_mm256_storeu_ps(c+3*stride+8,c31);
}
static void rows(const float*A,const float*B,float*E,int row0){
    memset(E,0,(size_t)MC*N*sizeof(float));
    for(int pc=0;pc<N;pc+=KC){
        for(int r=0;r<MC;r+=MR)for(int k=0;k<KC;++k)for(int i=0;i<MR;++i)
            apack[((r/MR)*KC+k)*MR+i]=A[(size_t)(row0+r+i)*N+pc+k];
        for(int jc=0;jc<N;jc+=NC){
            for(int q=0;q<NC;q+=NR)for(int k=0;k<KC;++k)for(int j=0;j<NR;++j)
                bpack[((q/NR)*KC+k)*NR+j]=B[(size_t)(pc+k)*N+jc+q+j];
            for(int r=0;r<MC;r+=MR)for(int q=0;q<NC;q+=NR)
                micro(apack+(size_t)r*KC,bpack+(size_t)q*KC,E+(size_t)r*N+jc+q,N);
        }
    }
}
}

int main(){
    constexpr int N=8192, ROWS=4096;
    float *A=nullptr,*B=nullptr,*E=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**)&E,4096,(size_t)128*N*4))return 2;
    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);
    volatile double checksum=0;
    for(int r=0;r<ROWS;r+=128){
        check_float::rows(A,B,E,r);
        for(int i=0;i<128;++i)checksum+=E[(size_t)i*N]+E[(size_t)i*N+N-1];
    }
    for(int i=0;i<4;++i){
        int row=ROWS-128+i*31,col=i*1701;
        double expected=0;
        for(int kk=0;kk<N;++kk)expected+=(double)A[(size_t)row*N+kk]*B[(size_t)kk*N+col];
        if(!isfinite(E[(size_t)(row-(ROWS-128))*N+col])||fabs(expected-E[(size_t)(row-(ROWS-128))*N+col])>0.1)return 3;
    }
    printf("checksum=%.12g\n",(double)checksum);
    free(E);free(B);free(A);return 0;
}

评测结果