【发布时间】：2014-07-10 13:08:22
【问题描述】：
问题
我正在研究 OpenBLAS 或 GotoBLAS 等高性能矩阵乘法算法，我正在尝试重现一些结果。这个问题涉及矩阵乘法算法的内核。具体来说，我正在研究计算C += AB，其中A 和B 是double 类型的2x2 矩阵，在我的CPU 的峰值速度下。有两种方法可以做到这一点。一种方法是使用 SIMD 指令。第二种方法是使用 SIMD 寄存器直接在汇编中编码。
到目前为止我所看到的
所有相关论文、课程网页、许多与该主题相关的 SO Q&A（太多无法列出），我在我的计算机上编译了 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 的峰值。作为参考，当两个内核都在峰值运行时，Intel Turbo Boost 提供了 2 步加速，我应该得到 22.4 DP Gflops 的理论峰值。
设置
我将我的 2x2 矩阵声明为 double 并使用随机条目初始化它们，如下面的代码 sn-p 所示。
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。代码 sn -p blow 说明了这部分。
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 个doubles。这是我在内核中进行 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);
}
汇编（英特尔语法）
我的第一次尝试是为此部件创建一个单独的汇编例程，并从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
问题
-
如果内联 ASM 和 SIMD 实现都产生相同的程序集输出，为什么程序集版本要慢得多？就好像内联程序集没有内联，第二组结果显示“内联”ASM 与“非内联”SIMD 的性能相同，这一点很明显。我能找到的唯一解释是在Agner Fog Volume 2 page 6：
编译代码可能比汇编代码更快，因为编译器可以使 程序间优化和整个程序优化。大会 程序员通常必须使用定义明确的调用来制作定义明确的函数 遵循所有调用约定的接口，以使代码可测试和 可验证的。这会阻止编译器使用的许多优化方法，例如 作为函数内联、寄存器分配、常量传播、公共子表达式 跨职能消除、跨职能调度等。这些 可以通过使用带有内部函数的 C++ 代码而不是 汇编代码。
但两个版本的汇编输出完全相同。
为什么我在第一组结果中看到 44 Gflops？这远高于我计算的 12 Gflops 峰值，如果我以单精度计算运行两个内核，这也是我所期望的。
编辑 1
评论说可能有死代码消除我可以确认 SIMd 指令正在发生这种情况。 -S 输出显示 SIMD 的 for 循环仅将 C 矩阵归零。我可以通过使用-O0 关闭编译器优化来禁用它。在这种情况下，SIMD 的运行速度是 ASM 的 3 倍，但 ASM 仍然以完全相同的速度运行。现在标准也是非零的，但在 10^-16 时仍然可以。我还看到内联 ASM 版本与 APP 和 NO_APP 标记内联，但它也在 for 循环中展开了 8 次。我认为展开多次会严重影响性能，因为我通常展开循环 4 次。根据我的经验，再多的东西似乎会降低性能。
【问题讨论】：
-
当你内联函数时，GCC 可能会跳过循环中的函数，因为 `C[0] = 0.0; C[1] = 0.0； C[2] = 0.0； C[3] = 0.0;` 您仍然可以在程序集中看到该函数，因为您在规范之前调用它。
-
使用 2x2 矩阵无法接近峰值性能。您需要每个时钟周期进行一次乘法、一次加法和一次加载。但是由于您还必须从 A 矩阵中读取数据，所以它更像是两次加载、一次乘法和一次加法。在我的代码中，我使用 AVX 执行 64x64 块。我做了 9 次加载（1 次来自 A，8 次来自 B）、8 次乘法和每行 8 次加法，因此我更接近每个时钟周期的一次乘法、一次加法和一次加载。 n 越大越好，除了 nxn 矩阵也需要适合 L1 缓存。
-
不知道你是否也有机会使用
FORTRAN。 -
我也有一个 4x4 版本，但在整个矩阵乘法框架中，它的运行速度比使用 2x2 时慢很多。不知道为什么...ja72 我不知道这有多大帮助，gcc 可能会将其编译为具有正确优化的相同机器代码。
-
@matmul，4x4 矩阵的失败次数是 2x2 的 8 倍，所以它应该花费更多时间，但也许你的意思是 GFLOPS 较低。那会很奇怪。您将读取 A 一次和 B 两次的每一行，然后执行 2 add 和 2 mult：3 加载，4 计算。其中 2x2 是 2 负载，2 计算，因此负载与计算的比率对于 2x2 更高，但我承认我还没有研究过这么小的矩阵的性能。