代码

#define BLOCK_SIZE_L1 32
#define BLOCK_SIZE_L2 64
#define BLOCK_SIZE_L3 128
#define MIN(a,b) ((a)<(b)?(a):(b))
#pragma GCC optimize("-Ofast","-funroll-all-loops","-ffast-math","-ftree-vectorize")
__attribute__((optimize("tree-vectorize"))) void matrix_multiply(int n_,const double* B,const double* A,double* C)
{
	const int n=1024;
	__attribute__((aligned(64))) double c[BLOCK_SIZE_L1][BLOCK_SIZE_L3];
	for(int ii=0;ii<n;ii+=BLOCK_SIZE_L3)
		for(int kk=0;kk<n;kk+=BLOCK_SIZE_L1)
		{
			#define ui BLOCK_SIZE_L3
			#define uj BLOCK_SIZE_L2
			#define uk BLOCK_SIZE_L1
			for(int k=0;k<uk;++k)
				for(int i=0;i<ui;++i) c[k][i]=0;
			for(int jj=0;jj<n;jj+=BLOCK_SIZE_L2)
			{
				const double*__restrict__ bb=B+jj+kk*n;
				const double*__restrict__ aa=A+ii+jj*n;
				bb=(const double*)((long long)bb&(~4096));
				aa=(const double*)((long long)aa&(~4096));
				for(int k=0;k<uk;++k)
					for(int j=0;j<uj;j++)
						for(int i=0;i<ui;i++)
							c[k][i]+=aa[i+j*n]*bb[j+k*n];
			}
			double*__restrict__ cc=C+ii+kk*n;
			cc=(double*)((long long)cc&(~4096));
			for(int k=0;k<uk;++k)
				for(int i=0;i<ui;++i)
					cc[i+k*n]+=c[k][i];
		}
}

评测结果