最多不能超过50％。 矩阵乘法的理论性能
问题
我正在学习HPC和代码优化。 我试图在Goto的开创性矩阵乘法论文（ http://www.cs.utexas.edu/users/pingali/CS378/2008sp/papers/gotoPaper.pdf ）中复制结果。 尽管我付出了最大的努力，但我无法超过理论CPU最高性能的50％。
背景
请参阅此处的相关问题（ 优化的2×2矩阵乘法：慢速assembly与快速SIMD ），包括有关我的硬件的信息
我尝试过的
这篇相关论文（ http://www.cs.utexas.edu/users/flame/pubs/blis3_ipdps14.pdf ）对Goto的算法结构有很好的描述。 我在下面提供了我的源代码。
我的问题
我要求一般帮助。 我一直在研究这个问题太久了，尝试了很多不同的算法，内联汇编，各种尺寸的内核（2×2,4×4,2×8，…，mxn m和n大），但我似乎无法打破50％CPU Gflops 。 这纯粹是出于教育目的，而不是作业。
源代码
希望是可以理解的。 请问是否。 我按照上面第2篇文章中的描述设置了宏结构（for loops）。 我按照两篇论文中的讨论打包矩阵，并在图11中以图形方式显示（ http://www.cs.utexas.edu/users/flame/pubs/BLISTOMSrev2.pdf ）。 我的内核计算2×8块，因为这似乎是Nehalem架构的最佳计算（参见GotoBLAS源代码 – 内核）。 内核基于计算排名1更新的概念，如此处所述（ http://code.google.com/p/blis/source/browse/config/template/kernels/3/bli_gemm_opt_mxn.c ）
#include #include #include #include #include #include #include #include // define some prefetch functions #define PREFETCHNTA(addr,nrOfBytesAhead) \ _mm_prefetch(((char *)(addr))+nrOfBytesAhead,_MM_HINT_NTA) #define PREFETCHT0(addr,nrOfBytesAhead) \ _mm_prefetch(((char *)(addr))+nrOfBytesAhead,_MM_HINT_T0) #define PREFETCHT1(addr,nrOfBytesAhead) \ _mm_prefetch(((char *)(addr))+nrOfBytesAhead,_MM_HINT_T1) #define PREFETCHT2(addr,nrOfBytesAhead) \ _mm_prefetch(((char *)(addr))+nrOfBytesAhead,_MM_HINT_T2) // define a min function #ifndef min #define min( a, b ) ( ((a) < (b)) ? (a) : (b) ) #endif // zero a matrix void zeromat(double *C, int n) { int i = n; while (i--) { int j = n; while (j--) { *(C + i*n + j) = 0.0; } } } // compute a 2x8 block from (2 x kc) x (kc x 8) matrices inline void __attribute__ ((gnu_inline)) __attribute__ ((aligned(64))) dgemm_2x8_sse( int k, const double* restrict a1, const int cs_a, const double* restrict b1, const int rs_b, double* restrict c11, const int rs_c ) { register __m128d xmm1, xmm4, // r8, r9, r10, r11, r12, r13, r14, r15; // accumulators // 10 registers declared here r8 = _mm_xor_pd(r8,r8); // ab r9 = _mm_xor_pd(r9,r9); r10 = _mm_xor_pd(r10,r10); r11 = _mm_xor_pd(r11,r11); r12 = _mm_xor_pd(r12,r12); // ab + 8 r13 = _mm_xor_pd(r13,r13); r14 = _mm_xor_pd(r14,r14); r15 = _mm_xor_pd(r15,r15); // PREFETCHT2(b1,0); // PREFETCHT2(b1,64); //int l = k; while (k--) { //PREFETCHT0(a1,0); // fetch 64 bytes from a1 // i = 0 xmm1 = _mm_load1_pd(a1); xmm4 = _mm_load_pd(b1); xmm4 = _mm_mul_pd(xmm1,xmm4); r8 = _mm_add_pd(r8,xmm4); xmm4 = _mm_load_pd(b1 + 2); xmm4 = _mm_mul_pd(xmm1,xmm4); r9 = _mm_add_pd(r9,xmm4); xmm4 = _mm_load_pd(b1 + 4); xmm4 = _mm_mul_pd(xmm1,xmm4); r10 = _mm_add_pd(r10,xmm4); xmm4 = _mm_load_pd(b1 + 6); xmm4 = _mm_mul_pd(xmm1,xmm4); r11 = _mm_add_pd(r11,xmm4); // // i = 1 xmm1 = _mm_load1_pd(a1 + 1); xmm4 = _mm_load_pd(b1); xmm4 = _mm_mul_pd(xmm1,xmm4); r12 = _mm_add_pd(r12,xmm4); xmm4 = _mm_load_pd(b1 + 2); xmm4 = _mm_mul_pd(xmm1,xmm4); r13 = _mm_add_pd(r13,xmm4); xmm4 = _mm_load_pd(b1 + 4); xmm4 = _mm_mul_pd(xmm1,xmm4); r14 = _mm_add_pd(r14,xmm4); xmm4 = _mm_load_pd(b1 + 6); xmm4 = _mm_mul_pd(xmm1,xmm4); r15 = _mm_add_pd(r15,xmm4); a1 += cs_a; b1 += rs_b; //PREFETCHT2(b1,0); //PREFETCHT2(b1,64); } // copy result into C PREFETCHT0(c11,0); xmm1 = _mm_load_pd(c11); xmm1 = _mm_add_pd(xmm1,r8); _mm_store_pd(c11,xmm1); xmm1 = _mm_load_pd(c11 + 2); xmm1 = _mm_add_pd(xmm1,r9); _mm_store_pd(c11 + 2,xmm1); xmm1 = _mm_load_pd(c11 + 4); xmm1 = _mm_add_pd(xmm1,r10); _mm_store_pd(c11 + 4,xmm1); xmm1 = _mm_load_pd(c11 + 6); xmm1 = _mm_add_pd(xmm1,r11); _mm_store_pd(c11 + 6,xmm1); c11 += rs_c; PREFETCHT0(c11,0); xmm1 = _mm_load_pd(c11); xmm1 = _mm_add_pd(xmm1,r12); _mm_store_pd(c11,xmm1); xmm1 = _mm_load_pd(c11 + 2); xmm1 = _mm_add_pd(xmm1,r13); _mm_store_pd(c11 + 2,xmm1); xmm1 = _mm_load_pd(c11 + 4); xmm1 = _mm_add_pd(xmm1,r14); _mm_store_pd(c11 + 4,xmm1); xmm1 = _mm_load_pd(c11 + 6); xmm1 = _mm_add_pd(xmm1,r15); _mm_store_pd(c11 + 6,xmm1); } // packs a matrix into rows of slivers inline void __attribute__ ((gnu_inline)) __attribute__ ((aligned(64))) rpack( double* restrict dst, const double* restrict src, const int kc, const int mc, const int mr, const int n) { double tmp[mc*kc] __attribute__ ((aligned(64))); double* restrict ptr = &tmp[0]; for (int i = 0; i < mc; ++i) for (int j = 0; j < kc; ++j) *ptr++ = *(src + i*n + j); ptr = &tmp[0]; //const int inc_dst = mr*kc; for (int k = 0; k < mc; k+=mr) for (int j = 0; j < kc; ++j) for (int i = 0; i < mr*kc; i+=kc) *dst++ = *(ptr + k*kc + j + i); } // packs a matrix into columns of slivers inline void __attribute__ ((gnu_inline)) __attribute__ ((aligned(64))) cpack(double* restrict dst, const double* restrict src, const int nc, const int kc, const int nr, const int n) { double tmp[kc*nc] __attribute__ ((aligned(64))); double* restrict ptr = &tmp[0]; for (int i = 0; i < kc; ++i) for (int j = 0; j < nc; ++j) *ptr++ = *(src + i*n + j); ptr = &tmp[0]; // const int inc_k = nc/nr; for (int k = 0; k < nc; k+=nr) for (int j = 0; j < kc*nc; j+=nc) for (int i = 0; i < nr; ++i) *dst++ = *(ptr + k + i + j); } void blis_dgemm_ref( const int n, const double* restrict A, const double* restrict B, double* restrict C, const int mc, const int nc, const int kc ) { int mr = 2; int nr = 8; double locA[mc*kc] __attribute__ ((aligned(64))); double locB[kc*nc] __attribute__ ((aligned(64))); int ii,jj,kk,i,j; #pragma omp parallel num_threads(4) shared(A,B,C) private(ii,jj,kk,i,j,locA,locB) {//use all threads in parallel #pragma omp for // partitions C and B into wide column panels for ( jj = 0; jj < n; jj+=nc) { // A and the current column of B are partitioned into col and row panels for ( kk = 0; kk < n; kk+=kc) { cpack(locB, B + kk*n + jj, nc, kc, nr, n); // partition current panel of A into blocks for ( ii = 0; ii < n; ii+=mc) { rpack(locA, A + ii*n + kk, kc, mc, mr, n); for ( i = 0; i < min(n-ii,mc); i+=mr) { for ( j = 0; j < min(n-jj,nc); j+=nr) { // inner kernel that compues 2 x 8 block dgemm_2x8_sse( kc, locA + i*kc , mr, locB + j*kc , nr, C + (i+ii)*n + (j+jj), n ); } } } } } } } double compute_gflops(const double time, const int n) { // computes the gigaflops for a square matrix-matrix multiplication double gflops; gflops = (double) (2.0*n*n*n)/time/1.0e9; return(gflops); } // ******* MAIN ********// void main() { clock_t time1, time2; double time3; double gflops; const int trials = 10; int nmax = 4096; printf("%10s %10s\n","N","Gflops/s"); int mc = 128; int kc = 256; int nc = 128; for (int n = kc; n <= nmax; n+=kc) { //assuming kc is the max dim double *A = NULL; double *B = NULL; double *C = NULL; A = _mm_malloc (n*n * sizeof(*A),64); B = _mm_malloc (n*n * sizeof(*B),64); C = _mm_malloc (n*n * sizeof(*C),64); srand(time(NULL)); // Create the matrices for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { A[i*n + j] = (double) rand()/RAND_MAX; B[i*n + j] = (double) rand()/RAND_MAX; //D[j*n + i] = B[i*n + j]; // Transpose C[i*n + j] = 0.0; } } // warmup zeromat(C,n); blis_dgemm_ref(n,A,B,C,mc,nc,kc); zeromat(C,n); time2 = 0; for (int count = 0; count < trials; count++){// iterations per experiment here time1 = clock(); blis_dgemm_ref(n,A,B,C,mc,nc,kc); time2 += clock() - time1; zeromat(C,n); } time3 = (double)(time2)/CLOCKS_PER_SEC/trials; gflops = compute_gflops(time3, n); printf("%10d %10f\n",n,gflops); _mm_free(A); _mm_free(B); _mm_free(C); } printf("tests are done\n"); } 编辑1
OS = Win 7 64位
编译器= gcc 4.8.1，但32位和mingw（也是32位。我正在努力获得mingw64的“非安装”版本，因此我可以使用更多XMM寄存器生成更快的代码/工作，等等。如果有人有链接到一个mingw64安装类似于mingw-get请发布。我的工作计算机有太多的管理限制。
填料
您似乎经常打包A矩阵的块。 你做
rpack(locA, A + ii*n + kk, kc, mc, mr, n);
但这只取决于ii和kk而不是jj但它在jj的内部循环中，所以你为jj每次迭代重新包装相同的东西。 我不认为这是必要的。 在我的代码中，我在矩阵乘法之前进行打包。 可能在矩阵乘法中打包更有效，而值仍然在缓存中，但这样做比较棘手。 但是打包是一个O（n ^ 2）运算，矩阵乘法是一个O（n ^ 3）运算，因此在大矩阵的矩阵乘法之外打包并不是非常低效（我也知道从测试开始 – 注释掉包装只会将效率提高几个百分点）。 但是，通过使用rpack重新打包每个jj迭代，您已经有效地使其成为O（n ^ 3）操作。
壁垒时间
你想要墙上的时间。 在Unix上，clock（）函数不会返回挂起时间 （尽管它在带有MSVC的Windows上）。 它返回每个线程的累积时间。 这是我在OpenMP上看到的最常见错误之一。
使用omp_get_wtime()来获取omp_get_wtime()时间。
请注意，我不知道clock()函数如何与MinGW或MinGW-w64一起使用（它们是单独的项目）。 MinGW链接到MSVCRT，所以我猜想minGW的clock()返回与MSVC一样的挂起时间。 但是，MinGW-w64没有链接到MSVCRT（据我所知它链接到像glibc这样的东西）。 MinGW-w64中的clock()与使用Unix的clock()执行相同的操作。
超线程
超线程适用于经常使CPU停顿的代码。 这实际上是大多数代码，因为编写不会使CPU停滞的代码非常困难。 这就是英特尔发明超线程的原因。 任务切换更容易，并为CPU提供别的东西，而不是优化代码。 但是，对于高度优化的代码，超线程实际上可能会产生更糟糕的结果。 在我自己的矩阵乘法代码中，情况肯定如此。 将线程数设置为您拥有的物理核心数（在您的情况下为两个）。
我的代码
以下是我的代码。 我没有在这里包含inner64函数。 您可以在MSVC和GCC之间的性能差异中找到它， 以获得高度优化的矩阵乘法代码 （ AddDot4x4_vec_block_8wide的令人讨厌和误导性的名称）
我在阅读Goto论文之前以及在阅读Agner Fog的优化手册之前编写了这段代码。 您似乎在主循环中重新排序/打包矩阵。 这可能更有意义。 我不认为我按照你的方式重新排序它们，而且我只重新排序其中一个输入矩阵（B），而不是像你一样。
在Linux和GCC上，我的系统（Xeon E5-1620@3.6）上的此代码的性能约为此矩阵大小（4096×4096）峰值的75％。 对于这种矩阵大小，英特尔的MKL在我的系统中达到峰值的94％左右，因此显然还有改进的余地。
#include #include #include #include extern "C" void inner64(const float *a, const float *b, float *c); void (*fp)(const float *a, const float *b, float *c) = inner64; void reorder(float * __restrict a, float * __restrict b, int n, int bs) { int nb = n/bs; #pragma omp parallel for for(int i=0; i