代码

#pragma GCC optimize("Ofast")
#pragma GCC target("avx")
#pragma GCC target("avx2")
#pragma GCC target("fma")
#include"bits/stdc++.h"
#include<immintrin.h>
#include<avxintrin.h>
#include<avx2intrin.h>
#include<fmaintrin.h>
#define f(_b) sum0=__builtin_ia32_vfmaddpd256(*(a1+(_b)),*(b1+(_b)),sum0)
#define re8(_b) f(_b),f(_b+1),f(_b+2),f(_b+3),f(_b+4),f(_b+5),f(_b+6),f(_b+7)
#define re16(_b) re8(_b),re8(_b+8)
#define re64(_b) re16(_b),re16(_b+16),re16(_b+32),re16(_b+48)
#define re256(_b) re64(_b),re64(_b+64),re64(_b+128),re64(_b+192)
#define re1024(_b) re256(_b),re256(_b+256),re256(_b+512),re256(_b+768)
#define re4096(_b) re1024(_b),re1024(_b+1024),re1024(_b+2048),re1024(_b+3072)
#define re8192(_b) re4096(_b),re4096(_b+4096)
#define N 1024
using namespace std;
double tmp[1048576];
void matrix_multiply(int n, const double *a, const double *b, double *c)
{
	for(int i=0;i<N;i++)for(int j=0;j<N;j++)*(tmp+(i<<10)+j)=*(b+(j<<10)+i);
	for(int i=0;i<N;i++)for(int j=0;j<N;j++)
	{
		__v4df sum0((__v4df){0,0,0,0}),*a1=(__v4df*)(a+(i<<10)),*b1=(__v4df*)(tmp+(j<<10));
		re256(0);
		*(c++)=sum0[0]+sum0[1]+sum0[2]+sum0[3];
	}
}

评测结果