代码

#define BLOCK_SIZE_L1 32
#define BLOCK_SIZE_L2 64
#define BLOCK_SIZE_L3 128
void square_dgemm(int n,const double* A,const double* B,double* C)
{
	__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)
		{
			int ui=n-ii,uk=n-kk;
			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)
			{
				double*__restrict__ bb=B+jj+kk*n;
				double*__restrict__ aa=A+ii+jj*n;
				int uj=n-jj;
				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;
			for(int k=0;k<uk;++k)
				for(int i=0;i<ui;++i)
					cc[i+k*n]+=c[k][i];
		}
}

评测结果