如何使用Dot产品获得峰值CPU性能？
问题
我一直在研究HPC，特别是使用矩阵乘法作为我的项目（参见我在配置文件中的其他post）。 我在那些方面取得了不错的成绩，但还不够好。 我退后一步看看我能用点积计算做得多好。
点乘积与矩阵乘法
点积更简单，我可以在不处理打包和其他相关问题的情况下测试HPC概念。 缓存阻塞仍然是一个问题，这是我的第二个问题。
算法
将两个double数组A和B n相应元素相乘并求它们。 assembly中的double点产品只是一系列movapd ， mulpd ， addpd 。 以巧妙的方式展开和排列，可以使movapd / mulpd / addpd组在不同的xmm寄存器上运行，因此是独立的，优化了流水线操作。 当然，事实certificate，这与我的CPU无序执行无关。 还要注意，重新安排需要剥离最后一次迭代。
其他假设
我不是在编写通用点积的代码。 代码是针对特定尺寸的，我不处理边缘情况。 这只是为了测试HPC概念并查看我可以获得的CPU使用类型。
结果
编译为gcc -std=c99 -O2 -m32 -mincoming-stack-boundary=2 -msse3 -mfpmath=sse,387 -masm=intel 。 我和平时不同的电脑。 这台计算机有一个i5 540m ，在两步Intel Turbo Boost之后， 2.8 GHz * 4 FLOPS/cycle/core = 11.2 GFLOPS/s per core可以获得2.8 GHz * 4 FLOPS/cycle/core = 11.2 GFLOPS/s per core （这两个核心现在都处于使用状态，所以它只有2步… a如果我关闭一个核心，可以进行4步增强。 当设置为使用一个线程运行时，32位LINPACK大约为9.5 GFLOPS / s。
N Total Gflops/s Residual 256 5.580521 1.421085e-014 384 5.734344 -2.842171e-014 512 5.791168 0.000000e+000 640 5.821629 0.000000e+000 768 5.814255 2.842171e-014 896 5.807132 0.000000e+000 1024 5.817208 -1.421085e-013 1152 5.805388 0.000000e+000 1280 5.830746 -5.684342e-014 1408 5.881937 -5.684342e-014 1536 5.872159 -1.705303e-013 1664 5.881536 5.684342e-014 1792 5.906261 -2.842171e-013 1920 5.477966 2.273737e-013 2048 5.620931 0.000000e+000 2176 3.998713 1.136868e-013 2304 3.370095 -3.410605e-013 2432 3.371386 -3.410605e-013
问题1
我怎么能比这更好？ 我甚至没有接近峰值表现。 我已经将assembly代码优化到了高天堂。 进一步展开可能会增加一点，但较少的展开似乎会降低性能。
问题2
当n > 2048 ，您可以看到性能下降。 这是因为我的L1缓存是32KB，当n = 2048且A和B是double ，它们占用整个缓存。 任何更大的，他们从记忆流。
我试过缓存阻塞（未在源代码中显示），但也许我做错了。 任何人都可以提供一些代码或解释如何阻止缓存的点积？
源代码
#include #include #include #include #include #include #include #include #include // computes 8 dot products #define KERNEL(address) \ "movapd xmm4, XMMWORD PTR [eax+"#address"] \n\t" \ "mulpd xmm7, XMMWORD PTR [edx+48+"#address"] \n\t" \ "addpd xmm2, xmm6 \n\t" \ "movapd xmm5, XMMWORD PTR [eax+16+"#address"] \n\t" \ "mulpd xmm4, XMMWORD PTR [edx+"#address"] \n\t" \ "addpd xmm3, xmm7 \n\t" \ "movapd xmm6, XMMWORD PTR [eax+96+"#address"] \n\t" \ "mulpd xmm5, XMMWORD PTR [edx+16+"#address"] \n\t" \ "addpd xmm0, xmm4 \n\t" \ "movapd xmm7, XMMWORD PTR [eax+112+"#address"] \n\t" \ "mulpd xmm6, XMMWORD PTR [edx+96+"#address"] \n\t" \ "addpd xmm1, xmm5 \n\t" #define PEELED(address) \ "movapd xmm4, XMMWORD PTR [eax+"#address"] \n\t" \ "mulpd xmm7, [edx+48+"#address"] \n\t" \ "addpd xmm2, xmm6 \n\t" \ "movapd xmm5, XMMWORD PTR [eax+16+"#address"] \n\t" \ "mulpd xmm4, XMMWORD PTR [edx+"#address"] \n\t" \ "addpd xmm3, xmm7 \n\t" \ "mulpd xmm5, XMMWORD PTR [edx+16+"#address"] \n\t" \ "addpd xmm0, xmm4 \n\t" \ "addpd xmm1, xmm5 \n\t" inline double __attribute__ ((gnu_inline)) __attribute__ ((aligned(64))) ddot_ref( int n, const double* restrict A, const double* restrict B) { double sum0 = 0.0; double sum1 = 0.0; double sum2 = 0.0; double sum3 = 0.0; double sum; for(int i = 0; i < n; i+=4) { sum0 += *(A + i ) * *(B + i ); sum1 += *(A + i+1) * *(B + i+1); sum2 += *(A + i+2) * *(B + i+2); sum3 += *(A + i+3) * *(B + i+3); } sum = sum0+sum1+sum2+sum3; return(sum); } inline double __attribute__ ((gnu_inline)) __attribute__ ((aligned(64))) ddot_asm ( int n, const double* restrict A, const double* restrict B) { double sum; __asm__ __volatile__ ( "mov eax, %[A] \n\t" "mov edx, %[B] \n\t" "mov ecx, %[n] \n\t" "pxor xmm0, xmm0 \n\t" "pxor xmm1, xmm1 \n\t" "pxor xmm2, xmm2 \n\t" "pxor xmm3, xmm3 \n\t" "movapd xmm6, XMMWORD PTR [eax+32] \n\t" "movapd xmm7, XMMWORD PTR [eax+48] \n\t" "mulpd xmm6, XMMWORD PTR [edx+32] \n\t" "sar ecx, 7 \n\t" "sub ecx, 1 \n\t" // peel "L%=: \n\t" KERNEL(64 * 0) KERNEL(64 * 1) KERNEL(64 * 2) KERNEL(64 * 3) KERNEL(64 * 4) KERNEL(64 * 5) KERNEL(64 * 6) KERNEL(64 * 7) KERNEL(64 * 8) KERNEL(64 * 9) KERNEL(64 * 10) KERNEL(64 * 11) KERNEL(64 * 12) KERNEL(64 * 13) KERNEL(64 * 14) KERNEL(64 * 15) "lea eax, [eax+1024] \n\t" "lea edx, [edx+1024] \n\t" " \n\t" "dec ecx \n\t" "jnz L%= \n\t" // end loop " \n\t" KERNEL(64 * 0) KERNEL(64 * 1) KERNEL(64 * 2) KERNEL(64 * 3) KERNEL(64 * 4) KERNEL(64 * 5) KERNEL(64 * 6) KERNEL(64 * 7) KERNEL(64 * 8) KERNEL(64 * 9) KERNEL(64 * 10) KERNEL(64 * 11) KERNEL(64 * 12) KERNEL(64 * 13) KERNEL(64 * 14) PEELED(64 * 15) " \n\t" "addpd xmm0, xmm1 \n\t" // summing result "addpd xmm2, xmm3 \n\t" "addpd xmm0, xmm2 \n\t" // cascading add "movapd xmm1, xmm0 \n\t" // copy xmm0 "shufpd xmm1, xmm0, 0x03 \n\t" // shuffle "addsd xmm0, xmm1 \n\t" // add low qword "movsd %[sum], xmm0 \n\t" // mov low qw to sum : // outputs [sum] "=m" (sum) : // inputs [A] "m" (A), [B] "m" (B), [n] "m" (n) : //register clobber "memory", "eax","ecx","edx","edi", "xmm0","xmm1","xmm2","xmm3","xmm4","xmm5","xmm6","xmm7" ); return(sum); } int main() { // timers LARGE_INTEGER frequency, time1, time2; double time3; QueryPerformanceFrequency(&frequency); // clock_t time1, time2; double gflops; int nmax = 4096; int trials = 1e7; double sum, residual; FILE *f = fopen("soddot.txt","w+"); printf("%16s %16s %16s\n","N","Total Gflops/s","Residual"); fprintf(f,"%16s %16s %16s\n","N","Total Gflops/s","Residual"); for(int n = 256; n <= nmax; n += 128 ) { double* A = NULL; double* B = NULL; A = _mm_malloc(n*sizeof(*A), 64); if (!A) {printf("A failed\n"); return(1);} B = _mm_malloc(n*sizeof(*B), 64); if (!B) {printf("B failed\n"); return(1);} srand(time(NULL)); // create arrays for(int i = 0; i < n; ++i) { *(A + i) = (double) rand()/RAND_MAX; *(B + i) = (double) rand()/RAND_MAX; } // warmup sum = ddot_asm(n,A,B); QueryPerformanceCounter(&time1); // time1 = clock(); for (int count = 0; count < trials; count++){ // sum = ddot_ref(n,A,B); sum = ddot_asm(n,A,B); } QueryPerformanceCounter(&time2); time3 = (double)(time2.QuadPart - time1.QuadPart) / frequency.QuadPart; // time3 = (double) (clock() - time1)/CLOCKS_PER_SEC; gflops = (double) (2.0*n*trials)/time3/1.0e9; residual = ddot_ref(n,A,B) - sum; printf("%16d %16f %16e\n",n,gflops,residual); fprintf(f,"%16d %16f %16e\n",n,gflops,residual); _mm_free(A); _mm_free(B); } fclose(f); return(0); // successful completion }
编辑：assembly说明
点积只是两个数的sum += a[i]*b[i]的重复和： sum += a[i]*b[i] 。 sum必须在第一次迭代之前初始化为0 。 矢量化，你可以一次做2个和，必须在末尾求和： [sum0 sum1] = [a[i] a[i+1]]*[b[i] b[i+1]] ， sum = sum0 + sum1 。 在（英特尔）程序集中，这是3个步骤（初始化后）：
pxor xmm0, xmm0 // accumulator [sum0 sum1] = [0 0] movapd xmm1, XMMWORD PTR [eax] // load [a[i] a[i+1]] into xmm1 mulpd xmm1, XMMWORD PTR [edx] // xmm1 = xmm1 * [b[i] b[i+1]] addpd xmm0, xmm1 // xmm0 = xmm0 + xmm1
此时你没什么特别的，编译器可以想出这个。 通过展开代码足以使用所有可用的xmm寄存器（32位模式下的8个寄存器），通常可以获得更好的性能。 因此，如果您将其展开4次，则可以使用所有8个寄存器xmm0到xmm7 。 您将有4个累加器和4个寄存器用于存储movapd和addpd的结果。 同样，编译器可以提出这个。 真正的思考部分是尝试提出一种管道代码的方法，即使MOV / MUL / ADD组中的每条指令在不同的寄存器上运行，以便所有3条指令同时执行（通常情况下是大多数CPU）。 这就是你打败编译器的方式。 因此，您必须对4x展开的代码进行模式化，这可能需要提前加载向量并剥离第一次或最后一次迭代。 这就是KERNEL(address) 。 为方便起见，我制作了4x展开的流水线代码的宏。 这样我只需更改address就可以轻松地将其展开为4的倍数。 每个KERNEL计算8个点产品。
要回答您的整体问题，您无法使用点积来达到最佳性能。
问题是您的CPU可以在每个时钟周期执行一次128位加载，并且为了完成点积，每个时钟周期需要两个128位加载。
但它比大n更糟糕。 你的第二个问题的答案是点积是内存绑定而不是计算绑定，所以它不能与具有快速核的大n并行化。 这里更好地解释了为什么 – 矢量化 – 循环 – 没有 – 性能改进 。 这是与快速核心并行化的一个大问题。 我花了一段时间才弄明白这一点，但学习起来非常重要。
实际上很少有基本算法可以从快速核心上的并行化中充分受益。 就BLAS算法而言，只有Level-3算法（O（n ^ 3））如矩阵乘法才真正受益于并行化。 在慢速核心上情况更好，例如GPU和Xeon Phi，因为内存速度和核心速度之间的差异要小得多。
如果你想找到一个算法，它可以接近峰值触发器，用于小n尝试例如标量*向量或标量*向量的总和。 第一种情况应该在每个时钟周期进行一次加载，一次多次和一次存储，第二种情况应该在每个时钟周期进行一次多次，一次加载和一次加载。
我在Knoppix 7.3 32位的Core 2 Duo P9600@2.67GHz上测试了以下代码。 我得到标量积的峰值的75％和标量积之和的峰值的75％。 标量积的翻转/周期为2，标量乘积的总和为4。
用g++ -msse2 -O3 -fopenmp foo.cpp -ffast-math编译g++ -msse2 -O3 -fopenmp foo.cpp -ffast-math
#include #include #include #include void scalar_product(double * __restrict a, int n) { a = (double*)__builtin_assume_aligned (a, 64); double k = 3.14159; for(int i=0; i