优化的2×2矩阵乘法：慢速assembly与快速SIMD

问题

我正在研究高性能矩阵乘法算法，如OpenBLAS或GotoBLAS，我正在尝试重现一些结果。 这个问题涉及矩阵乘法算法的内核。 具体来说，我正在研究计算C += AB ，其中A和B是在我的CPU的峰值速度下类型为double 2×2矩阵。 有两种方法可以做到这一点。 一种方法是使用SIMD指令。 第二种方法是使用SIMD寄存器直接在汇编代码中编码。

到目前为止我看过的是什么

所有相关的论文，课程网页，许多SO Q＆As处理主题（太多无法列出），我在我的计算机上编译了OpenBLAS，查看了OpenBLAS，GotoBLAS和BLIS源代码，Agner的手册。

硬件

我的CPU是Intel i5 – 540M。 您可以在cpu-world.com上找到相关的CPUID信息。 微体系结构是Nehalem（westmere），因此理论上每循环每个核心可以计算4个双精度触发器。 我将只使用一个核心（没有OpenMP），因此，对于超线程关闭和4步Intel Turbo Boost，我应该看到( 2.533 Ghz + 4*0.133 Ghz ) * ( 4 DP flops/core/cycle ) * ( 1 core ) = 12.27 DP Gflops的峰值( 2.533 Ghz + 4*0.133 Ghz ) * ( 4 DP flops/core/cycle ) * ( 1 core ) = 12.27 DP Gflops 。 作为参考，两个内核都运行在峰值时，Intel Turbo Boost提供了两步加速，我应该获得22.4 DP Gflops的理论峰值。

建立

我将我的2×2矩阵声明为double并使用随机条目对其进行初始化，如下面的代码片段所示。

 srand(time(NULL)); const int n = 2; double A[n*n]; double B[n*n]; double C[n*n]; double T[n*n]; for(int i = 0; i < n*n; i++){ A[i] = (double) rand()/RAND_MAX; B[i] = (double) rand()/RAND_MAX; C[i] = 0.0; } 

我使用朴素矩阵矩阵乘法（如下所示）计算一个真实的答案，它允许我通过视觉或通过计算所有元素的L2范数来检查我的结果

 // "true" answer for(int i = 0; i < n; i++) for(int j = 0; j < n; j++) for(int k = 0; k < n; k++) T[i*n + j] += A[i*n + k]*B[k*n + j]; 

为了运行代码并获得Gflops的估计值，我将每个乘法函数调用一次以进行预热，然后在for循环中执行它for进行maxiter次数，确保每次在计算C += AB时将C矩阵归零。 for循环放在两个clock()语句中，用于估计Gflops。 代码片段说明了这一部分。

 C[0] = 0.0; C[1] = 0.0; C[2] = 0.0; C[3] = 0.0; mult2by2(A,B,C); //warmup time1 = clock(); for(int i = 0; i < maxiter; i++){ mult2by2(A,B,C); C[0] = 0.0; C[1] = 0.0; C[2] = 0.0; C[3] = 0.0; } time2 = clock() - time1; time3 = (double)(time2)/CLOCKS_PER_SEC; gflops = (double) (2.0*n*n*n)/time3/1.0e9*maxiter; mult2by2(A,B,C); // to compute the norm against T norm = L2norm(n,C,T); 

SIMD代码

我的CPU支持128位向量，所以我可以在每个向量中加入2个double s。 这是我在内核中进行2×2矩阵乘法的主要原因。 SIMD代码一次计算整行C

  inline void __attribute__ ((gnu_inline)) __attribute__ ((aligned(16))) mult2by2B( const double* restrict A, const double* restrict B, double* restrict C ) { register __m128d xmm0, xmm1, xmm2, xmm3, xmm4; xmm0 = _mm_load_pd(C); xmm1 = _mm_load1_pd(A); xmm2 = _mm_load_pd(B); xmm3 = _mm_load1_pd(A + 1); xmm4 = _mm_load_pd(B + 2); xmm1 = _mm_mul_pd(xmm1,xmm2); xmm2 = _mm_add_pd(xmm1,xmm0); xmm1 = _mm_mul_pd(xmm3,xmm4); xmm2 = _mm_add_pd(xmm1,xmm2); _mm_store_pd(C,xmm2); xmm0 = _mm_load_pd(C + 2); xmm1 = _mm_load1_pd(A + 2); xmm2 = _mm_load_pd(B); xmm3 = _mm_load1_pd(A + 3); //xmm4 = _mm_load_pd(B + 2); xmm1 = _mm_mul_pd(xmm1,xmm2); xmm2 = _mm_add_pd(xmm1,xmm0); xmm1 = _mm_mul_pd(xmm3,xmm4); xmm2 = _mm_add_pd(xmm1,xmm2); _mm_store_pd(C + 2,xmm2); } 

Assmebly（英特尔语法）

我的第一次尝试是为这个部分创建一个单独的汇编程序，并从main调用它。 但是，它非常慢，因为我无法内联extern函数。 我将程序集编写为内联汇编，如下所示。 它与 gcc -S -std=c99 -O3 -msse3 -ffast-math -march=nocona -mtune=nocona -funroll-all-loops -fomit-frame-pointer -masm=intel 。 根据我对Nehalem微体系结构图的理解，该处理器可以并行执行SSE ADD ， SSE MUL和SSE MOV ，这解释了MUL ， ADD ， MOV指令的交错。 您会注意到上面的SIMD说明顺序不同，因为我对Agner Fog的手册有不同的理解。 不过， gcc很聪明，上面的SIMD代码编译成内联版本中显示的程序集。

 inline void __attribute__ ((gnu_inline)) __attribute__ ((aligned(16))) mult2by2A ( const double* restrict A, const double* restrict B, double* restrict C ) { __asm__ __volatile__ ( "mov edx, %[A] \n\t" "mov ecx, %[B] \n\t" "mov eax, %[C] \n\t" "movapd xmm3, XMMWORD PTR [ecx] \n\t" "movapd xmm2, XMMWORD PTR [ecx+16] \n\t" "movddup xmm1, QWORD PTR [edx] \n\t" "mulpd xmm1, xmm3 \n\t" "addpd xmm1, XMMWORD PTR [eax] \n\t" "movddup xmm0, QWORD PTR [edx+8] \n\t" "mulpd xmm0, xmm2 \n\t" "addpd xmm0, xmm1 \n\t" "movapd XMMWORD PTR [eax], xmm0 \n\t" "movddup xmm4, QWORD PTR [edx+16] \n\t" "mulpd xmm4, xmm3 \n\t" "addpd xmm4, XMMWORD PTR [eax+16] \n\t" "movddup xmm5, QWORD PTR [edx+24] \n\t" "mulpd xmm5, xmm2 \n\t" "addpd xmm5, xmm4 \n\t" "movapd XMMWORD PTR [eax+16], xmm5 \n\t" : // no outputs : // inputs [A] "m" (A), [B] "m" (B), [C] "m" (C) : //register clobber "memory", "edx","ecx","eax", "xmm0","xmm1","xmm2","xmm3","xmm4","xmm5" ); } 

结果

我使用以下标志编译我的代码：

 gcc -std=c99 -O3 -msse3 -ffast-math -march=nocona -mtune=nocona -funroll-all-loops -fomit-frame-pointer -masm=intel 

maxiter = 1000000000的结果如下：

 ********** Inline ASM L2 norm: 0.000000e+000, Avg. CPU time: 9.563000, Avg. Gflops: 1.673115 ********** SIMD Version L2 norm: 0.000000e+000, Avg. CPU time: 0.359000, Avg. Gflops: 44.568245 

如果我强制SIMD版本没有内联__attribute__ ((noinline)) ，结果是：

 ********** Inline ASM L2 norm: 0.000000e+000, Avg. CPU time: 11.155000, Avg. Gflops: 1.434334 ********** SIMD Version L2 norm: 0.000000e+000, Avg. CPU time: 11.264000, Avg. Gflops: 1.420455 

问题

1. 如果内联ASM和SIMD实现产生相同的汇编输出，为什么汇编版本要慢得多？ 就像内联汇编没有内联一样，第二组结果显示“内联”ASM与“无内线”SIMD的性能相同。 我能找到的唯一解释是在Agner Fog第2卷第6页 ：

   编译代码可能比汇编代码更快，因为编译器可以进行过程间优化和整个程序优化。 汇编程序员通常必须使用定义良好的调用接口来定义良好定义的函数，该调用接口遵循所有调用约定，以使代码可测试和可validation。 这可以防止编译器使用的许多优化方法，例如函数内联，寄存器分配，常量传播，跨函数的公共子表达式消除，跨函数调度等。这些优点可以通过使用具有内部函数的C ++代码而不是汇编代码来获得。 。

   但两个版本的汇编输出完全相同。

2. 为什么我在第一组结果中看到44 Gflops？ 这高于我计算的12 Gflops峰值，如果我使用单精度计算运行两个核心，这就是我所期望的。

编辑1评论说可能有死代码消除我可以确认SIMd指令正在发生这种情况。 -S输出显示SIMD的for循环仅为零矩阵。 我可以通过使用-O0关闭编译器优化来禁用它。 在这种情况下，SIMD的运行速度是ASM的3倍，但ASM仍以完全相同的速度运行。 规范现在也非零，但它仍然可以在10 ^ -16。 我还看到内联ASM版本正在内联APP和NO_APP标签，但它也在for循环中展开了8次。 我认为多次展开会严重影响性能，因为我通常会循环4次循环。 根据我的经验，更多的东西似乎会降低性能。

Answer 1

由于线路，GCC正在使用内在函数mult2by2B优化你的内联函数

 C[0] = 0.0; C[1] = 0.0; C[2] = 0.0; C[3] = 0.0; 

如果没有这条线，那么从Coliru开始计算机需要2.9秒http://coliru.stacked-crooked.com/a/992304f5f672e257

这条线只需要0.000001 http://coliru.stacked-crooked.com/a/9722c39bb6b8590a

您也可以在程序集中看到这一点。 如果您将以下代码放入http://gcc.godbolt.org/，您将看到使用该行代码完全跳过该函数。

但是，当你内联程序集时，GCC不会优化函数mult2by2A ，即使它内联它（即使它内联它）。 您也可以在程序集中看到这一点。

 #include  #include  // SSE2 #include  inline void __attribute__ ((gnu_inline)) __attribute__ ((aligned(16))) mult2by2B( const double* __restrict A, const double* __restrict B, double* __restrict C ) { register __m128d xmm0, xmm1, xmm2, xmm3, xmm4; xmm0 = _mm_load_pd(C); xmm1 = _mm_load1_pd(A); xmm2 = _mm_load_pd(B); xmm3 = _mm_load1_pd(A + 1); xmm4 = _mm_load_pd(B + 2); xmm1 = _mm_mul_pd(xmm1,xmm2); xmm2 = _mm_add_pd(xmm1,xmm0); xmm1 = _mm_mul_pd(xmm3,xmm4); xmm2 = _mm_add_pd(xmm1,xmm2); _mm_store_pd(C,xmm2); xmm0 = _mm_load_pd(C + 2); xmm1 = _mm_load1_pd(A + 2); xmm2 = _mm_load_pd(B); xmm3 = _mm_load1_pd(A + 3); //xmm4 = _mm_load_pd(B + 2); xmm1 = _mm_mul_pd(xmm1,xmm2); xmm2 = _mm_add_pd(xmm1,xmm0); xmm1 = _mm_mul_pd(xmm3,xmm4); xmm2 = _mm_add_pd(xmm1,xmm2); _mm_store_pd(C + 2,xmm2); } int main() { double A[4], B[4], C[4]; int maxiter = 10000000; //int maxiter = 1000000000; double dtime; dtime = omp_get_wtime(); for(int i = 0; i < maxiter; i++){ mult2by2B(A,B,C); C[0] = 0.0; C[1] = 0.0; C[2] = 0.0; C[3] = 0.0; } dtime = omp_get_wtime() - dtime; printf("%f %f %f %f\n", C[0], C[1], C[2], C[3]); //gflops = (double) (2.0*n*n*n)/time3/1.0e9*maxiter; printf("time %f\n", dtime); }