问题
我正在研究高性能矩阵乘法算法，例如 OpenBLAS 或 GotoBLAS，我正在尝试重现一些结果。这个问题涉及矩阵乘法算法的内部内核。具体来说，我正在研究计算 C += AB ，其中 A和 B是 double 类型的 2x2 矩阵以我的 CPU 的峰值速度。有两种方法可以做到这一点。一种方法是使用 SIMD 指令。第二种方法是使用 SIMD 寄存器直接在汇编中编码。
到目前为止我看过什么
所有相关论文、类(class)网页、许多关于该主题的 SO Q&A(太多无法一一列举)，我在我的电脑上编译了 OpenBLAS，浏览了 OpenBLAS、GotoBLAS 和 BLIS 源代码，Agner 的手册。
硬件
我的 CPU 是 Intel i5 - 540M。您可以在 cpu-world.com 上找到相关的 CPUID 信息。微架构是 Nehalem (westmere)，因此理论上每个内核每个周期可以计算 4 个 double 触发器。我将只使用一个内核(不使用 OpenMP)，因此在关闭超线程和 4 步 Intel Turbo Boost 后，我应该会看到峰值 ( 2.533 Ghz + 4*0.133 Ghz ) * ( 4 DP flops/core/cycle ) * ( 1 core ) = 12.27 DP Gflops .作为引用，当两个内核都在峰值运行时，Intel Turbo Boost 提供了 2 步加速，我应该得到 22.4 DP Gflops 的理论峰值。 .
设置
我将我的 2x2 矩阵声明为 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 中执行它。循环 maxiter次，确保将 C 归零每次计算时的矩阵 C += AB . 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 位 vector ，所以我可以容纳 2
double每个 vector 中的 s。这是我在内核中进行 2x2 矩阵乘法的主要原因。 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);
}
组装(英特尔语法)
我的第一次尝试是为这部分创建一个单独的 assembly 例程，并从
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
问题
Compiled code may be faster than assembly code because compilers can make inter-procedural optimization and whole-program optimization. The assembly programmer usually has to make well-defined functions with a well-defined call interface that obeys all calling conventions in order to make the code testable and verifiable. This prevents many of the optimization methods that compilers use, such as function inlining, register allocation, constant propagation, common subexpression elimination across functions, scheduling across functions, etc. These advantages can be obtained by using C++ code with intrinsic functions instead of assembly code.
但是两个版本的汇编器输出完全相同。
编辑 1
评论说可能会消除死代码，我可以确认 SIMd 指令正在发生这种情况。
-S输出显示 for SIMD 循环仅零 C矩阵。我可以通过使用 -O0 关闭编译器优化来禁用它.在这种情况下，SIMD 的运行速度是 ASM 的 3 倍，但 ASM 仍然以完全相同的速度运行。范数现在也非零，但在 10^-16 处仍然可以。我还看到内联 ASM 版本正在与 APP 内联。和 NO_APP标签，但它也在 for 中展开了 8 次环形。我认为多次展开会严重影响性能，因为我通常会展开循环 4 次。根据我的经验，更多的东西似乎会降低性能。
最佳答案
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 , away (即使它内联它)。您也可以在程序集中看到这一点。#include <stdio.h>
#include <emmintrin.h> // SSE2
#include <omp.h>
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);
}
关于c - 优化的 2x2 矩阵乘法 : Slow assembly versus fast SIMD，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/23790842/