如何让 GCC 完全展开这个循环（即剥离这个循环）？
How to ask GCC to completely unroll this loop (i.e., peel this loop)?
有没有办法指示 GCC（我使用的是 4.8.4）完全展开底部函数中的 while 循环 ，即剥离这个环形？循环的迭代次数在编译时已知：58.
让我先解释一下我尝试过的方法。
通过检查 GAS 输出：
gcc -fpic -O2 -S GEPDOT.c
使用了 12 个寄存器 XMM0 - XMM11。如果我将标志 -funroll-loops 传递给 gcc:
gcc -fpic -O2 -funroll-loops -S GEPDOT.c
循环只展开了两次。我检查了 GCC 优化选项。 GCC 表示 -funroll-loops 也将打开 -frename-registers，因此当 GCC 展开循环时，其寄存器分配的优先选择是使用“剩余”寄存器。但是 XMM12 - XMM15 只剩下 4 个，所以 GCC 最多只能展开 2 次。如果有 48 个而不是 16 个可用的 XMM 寄存器，GCC 将毫无问题地展开 while 循环 4 次。
然而我又做了一个实验。我首先手动展开 while 循环两次，获得一个函数 GEPDOT_2。那么在
之间完全没有区别
gcc -fpic -O2 -S GEPDOT_2.c
和
gcc -fpic -O2 -funroll-loops -S GEPDOT_2.c
因为GEPDOT_2已经用完了所有寄存器，所以不执行展开。
GCC 确实会重命名寄存器以避免引入潜在的 虚假依赖。但我很确定我的 GEPDOT; 不会有这样的潜力。即使有，也不重要。我尝试自己展开循环，展开 4 次比展开 2 次快，比不展开快。当然我可以手动展开更多次，但是很繁琐。 GCC 可以为我做这个吗？谢谢。
// C file "GEPDOT.c"
#include <emmintrin.h>
void GEPDOT (double *A, double *B, double *C) {
__m128d A1_vec = _mm_load_pd(A); A += 2;
__m128d B_vec = _mm_load1_pd(B); B++;
__m128d C1_vec = A1_vec * B_vec;
__m128d A2_vec = _mm_load_pd(A); A += 2;
__m128d C2_vec = A2_vec * B_vec;
B_vec = _mm_load1_pd(B); B++;
__m128d C3_vec = A1_vec * B_vec;
__m128d C4_vec = A2_vec * B_vec;
B_vec = _mm_load1_pd(B); B++;
__m128d C5_vec = A1_vec * B_vec;
__m128d C6_vec = A2_vec * B_vec;
B_vec = _mm_load1_pd(B); B++;
__m128d C7_vec = A1_vec * B_vec;
A1_vec = _mm_load_pd(A); A += 2;
__m128d C8_vec = A2_vec * B_vec;
B_vec = _mm_load1_pd(B); B++;
int k = 58;
/* can compiler unroll the loop completely (i.e., peel this loop)? */
while (k--) {
C1_vec += A1_vec * B_vec;
A2_vec = _mm_load_pd(A); A += 2;
C2_vec += A2_vec * B_vec;
B_vec = _mm_load1_pd(B); B++;
C3_vec += A1_vec * B_vec;
C4_vec += A2_vec * B_vec;
B_vec = _mm_load1_pd(B); B++;
C5_vec += A1_vec * B_vec;
C6_vec += A2_vec * B_vec;
B_vec = _mm_load1_pd(B); B++;
C7_vec += A1_vec * B_vec;
A1_vec = _mm_load_pd(A); A += 2;
C8_vec += A2_vec * B_vec;
B_vec = _mm_load1_pd(B); B++;
}
C1_vec += A1_vec * B_vec;
A2_vec = _mm_load_pd(A);
C2_vec += A2_vec * B_vec;
B_vec = _mm_load1_pd(B); B++;
C3_vec += A1_vec * B_vec;
C4_vec += A2_vec * B_vec;
B_vec = _mm_load1_pd(B); B++;
C5_vec += A1_vec * B_vec;
C6_vec += A2_vec * B_vec;
B_vec = _mm_load1_pd(B);
C7_vec += A1_vec * B_vec;
C8_vec += A2_vec * B_vec;
/* [write-back] */
A1_vec = _mm_load_pd(C); C1_vec = A1_vec - C1_vec;
A2_vec = _mm_load_pd(C + 2); C2_vec = A2_vec - C2_vec;
A1_vec = _mm_load_pd(C + 4); C3_vec = A1_vec - C3_vec;
A2_vec = _mm_load_pd(C + 6); C4_vec = A2_vec - C4_vec;
A1_vec = _mm_load_pd(C + 8); C5_vec = A1_vec - C5_vec;
A2_vec = _mm_load_pd(C + 10); C6_vec = A2_vec - C6_vec;
A1_vec = _mm_load_pd(C + 12); C7_vec = A1_vec - C7_vec;
A2_vec = _mm_load_pd(C + 14); C8_vec = A2_vec - C8_vec;
_mm_store_pd(C,C1_vec); _mm_store_pd(C + 2,C2_vec);
_mm_store_pd(C + 4,C3_vec); _mm_store_pd(C + 6,C4_vec);
_mm_store_pd(C + 8,C5_vec); _mm_store_pd(C + 10,C6_vec);
_mm_store_pd(C + 12,C7_vec); _mm_store_pd(C + 14,C8_vec);
}
更新 1
感谢@user3386109 的评论，我想稍微扩展一下这个问题。 @user3386109 提出了一个 非常好的 问题。实际上，当有这么多并行指令要调度时，我确实怀疑编译器优化寄存器分配的能力。
我个人认为一个可靠的方法是先在 asm 内联汇编中编写循环体（这是 HPC 的关键），然后根据需要复制多次.今年早些时候我有一个不受欢迎的post：。代码略有不同，因为循环迭代次数 j 是一个函数参数，因此在编译时是未知的。在那种情况下我无法完全展开循环，所以我只复制了两次汇编代码，并将循环转换为标签并跳转。事实证明，我编写的程序集的最终性能比编译器生成的程序集高约 5%，这可能表明编译器未能以我们预期的最佳方式分配寄存器。
我曾经（现在仍然）是汇编编码的婴儿，所以这对我来说是一个很好的案例研究，可以让我学习一些 x86 汇编。但是在很长的运行中，我不倾向于编码GEPDOT 占很大比例的汇编。主要有以下三个原因：
- asm 内联汇编因为不可移植而受到批评。虽然我不明白为什么。也许是因为不同的机器有不同的寄存器被破坏？
- 编译器也越来越好。所以我还是更喜欢算法优化和更好的C编码习惯来帮助编译器产生好的输出；
- 最后一个原因更重要。迭代次数可能并不总是 58。我正在开发一个高性能矩阵分解子程序。对于缓存阻塞因子
nb，迭代次数将为 nb-2。我不会将 nb 作为函数参数，就像我在之前的 post 中所做的那样。这是一个机器特定参数，将被定义为一个宏。所以迭代次数在编译时是已知的，但可能因机器而异。猜猜我在为各种 nb 手动循环展开时需要做多少繁琐的工作。因此，如果有一种方法可以简单地指示编译器剥离一个循环，那就太好了。
如果您也能分享一些制作高性能但可移植的库的经验，我将不胜感激。
尝试调整优化器参数：
gcc -O3 -funroll-loops --param max-completely-peeled-insns=1000 --param max-completely-peel-times=100
这应该可以解决问题。
这不是答案，但其他尝试使用 GCC 向量化矩阵乘法的人可能会感兴趣。
下面，我假设 c 是行主顺序的 4×4 矩阵，a 是 4 行，n-列主矩阵（t运行sposed），b是一个4列， n-行矩阵行主序，运算为c = a × b + c，其中×表示矩阵乘法。
实现这个的最简单的函数是
void slow_4(double *c,
const double *a,
const double *b,
size_t n)
{
size_t row, col, i;
for (row = 0; row < 4; row++)
for (col = 0; col < 4; col++)
for (i = 0; i < n; i++)
c[4*row+col] += a[4*i+row] * b[4*i+col];
}
GCC 使用
为 SSE2/SSE3 生成了非常好的代码
#if defined(__SSE2__) || defined(__SSE3__)
typedef double vec2d __attribute__((vector_size (2 * sizeof (double))));
void fast_4(vec2d *c,
const vec2d *a,
const vec2d *b,
size_t n)
{
const vec2d *const b_end = b + 2L * n;
vec2d s00 = c[0];
vec2d s02 = c[1];
vec2d s10 = c[2];
vec2d s12 = c[3];
vec2d s20 = c[4];
vec2d s22 = c[5];
vec2d s30 = c[6];
vec2d s32 = c[7];
while (b < b_end) {
const vec2d b0 = b[0];
const vec2d b2 = b[1];
const vec2d a0 = { a[0][0], a[0][0] };
const vec2d a1 = { a[0][1], a[0][1] };
const vec2d a2 = { a[1][0], a[1][0] };
const vec2d a3 = { a[1][1], a[1][1] };
s00 += a0 * b0;
s10 += a1 * b0;
s20 += a2 * b0;
s30 += a3 * b0;
s02 += a0 * b2;
s12 += a1 * b2;
s22 += a2 * b2;
s32 += a3 * b2;
b += 2;
a += 2;
}
c[0] = s00;
c[1] = s02;
c[2] = s10;
c[3] = s12;
c[4] = s20;
c[5] = s22;
c[6] = s30;
c[7] = s32;
}
#endif
对于 AVX，GCC 可以做得更好
#if defined(__AVX__) || defined(__AVX2__)
typedef double vec4d __attribute__((vector_size (4 * sizeof (double))));
void fast_4(vec4d *c,
const vec4d *a,
const vec4d *b,
size_t n)
{
const vec4d *const b_end = b + n;
vec4d s0 = c[0];
vec4d s1 = c[1];
vec4d s2 = c[2];
vec4d s3 = c[3];
while (b < b_end) {
const vec4d bc = *(b++);
const vec4d ac = *(a++);
const vec4d a0 = { ac[0], ac[0], ac[0], ac[0] };
const vec4d a1 = { ac[1], ac[1], ac[1], ac[1] };
const vec4d a2 = { ac[2], ac[2], ac[2], ac[2] };
const vec4d a3 = { ac[3], ac[3], ac[3], ac[3] };
s0 += a0 * bc;
s1 += a1 * bc;
s2 += a2 * bc;
s3 += a3 * bc;
}
c[0] = s0;
c[1] = s1;
c[2] = s2;
c[3] = s3;
}
#endif
使用gcc-4.8.4(-O2 -march=x86-64 -mtune=generic -msse3)生成的程序集的SSE3版本本质上是
fast_4:
salq , %rcx
movapd (%rdi), %xmm13
addq %rdx, %rcx
cmpq %rcx, %rdx
movapd 16(%rdi), %xmm12
movapd 32(%rdi), %xmm11
movapd 48(%rdi), %xmm10
movapd 64(%rdi), %xmm9
movapd 80(%rdi), %xmm8
movapd 96(%rdi), %xmm7
movapd 112(%rdi), %xmm6
jnb .L2
.L3:
movddup (%rsi), %xmm5
addq , %rdx
movapd -32(%rdx), %xmm1
addq , %rsi
movddup -24(%rsi), %xmm4
movapd %xmm5, %xmm14
movddup -16(%rsi), %xmm3
movddup -8(%rsi), %xmm2
mulpd %xmm1, %xmm14
movapd -16(%rdx), %xmm0
cmpq %rdx, %rcx
mulpd %xmm0, %xmm5
addpd %xmm14, %xmm13
movapd %xmm4, %xmm14
mulpd %xmm0, %xmm4
addpd %xmm5, %xmm12
mulpd %xmm1, %xmm14
addpd %xmm4, %xmm10
addpd %xmm14, %xmm11
movapd %xmm3, %xmm14
mulpd %xmm0, %xmm3
mulpd %xmm1, %xmm14
mulpd %xmm2, %xmm0
addpd %xmm3, %xmm8
mulpd %xmm2, %xmm1
addpd %xmm14, %xmm9
addpd %xmm0, %xmm6
addpd %xmm1, %xmm7
ja .L3
.L2:
movapd %xmm13, (%rdi)
movapd %xmm12, 16(%rdi)
movapd %xmm11, 32(%rdi)
movapd %xmm10, 48(%rdi)
movapd %xmm9, 64(%rdi)
movapd %xmm8, 80(%rdi)
movapd %xmm7, 96(%rdi)
movapd %xmm6, 112(%rdi)
ret
生成的程序集(-O2 -march=x86-64 -mtune=generic -mavx)的AVX版本本质上是
fast_4:
salq , %rcx
vmovapd (%rdi), %ymm5
addq %rdx, %rcx
vmovapd 32(%rdi), %ymm4
cmpq %rcx, %rdx
vmovapd 64(%rdi), %ymm3
vmovapd 96(%rdi), %ymm2
jnb .L2
.L3:
addq , %rsi
vmovapd -32(%rsi), %ymm1
addq , %rdx
vmovapd -32(%rdx), %ymm0
cmpq %rdx, %rcx
vpermilpd [=14=], %ymm1, %ymm6
vperm2f128 [=14=], %ymm6, %ymm6, %ymm6
vmulpd %ymm0, %ymm6, %ymm6
vaddpd %ymm6, %ymm5, %ymm5
vpermilpd , %ymm1, %ymm6
vperm2f128 [=14=], %ymm6, %ymm6, %ymm6
vmulpd %ymm0, %ymm6, %ymm6
vaddpd %ymm6, %ymm4, %ymm4
vpermilpd [=14=], %ymm1, %ymm6
vpermilpd , %ymm1, %ymm1
vperm2f128 , %ymm6, %ymm6, %ymm6
vperm2f128 , %ymm1, %ymm1, %ymm1
vmulpd %ymm0, %ymm6, %ymm6
vmulpd %ymm0, %ymm1, %ymm0
vaddpd %ymm6, %ymm3, %ymm3
vaddpd %ymm0, %ymm2, %ymm2
ja .L3
.L2:
vmovapd %ymm5, (%rdi)
vmovapd %ymm4, 32(%rdi)
vmovapd %ymm3, 64(%rdi)
vmovapd %ymm2, 96(%rdi)
vzeroupper
ret
寄存器调度不是最优的，我猜，但它看起来也不糟糕。我个人对上面的内容很满意，此时没有尝试手动优化它。
在 Core i5-4200U 处理器（支持 AVX2）上，上述函数的快速版本在 1843 CPU 个 SSE3 周期（中位数）和 1248 中计算两个 4×256 矩阵的乘积AVX2 的周期。这归结为每个矩阵条目有 1.8 和 1.22 个周期。未向量化的慢速版本每个矩阵条目大约需要 11 个周期，以供比较。
（循环计数是中值，即一半的测试速度更快。我只是 运行 一些粗略的基准测试，大约 100k 次重复，所以请对这些数字持保留态度。）
在这个 CPU 上，缓存效果使得 AVX2 在 4×512 矩阵大小下仍为每个条目 1.2 个周期，但在 4×1024 时，它下降到 1.4，在 4×4096 时下降到 1.5 ，在 4×8192 到 1.8，以及在 4×65536 到 2.2 每个条目的周期。 SSE3 版本保持在每个条目 1.8 个周期，最高为 4×3072，此时它开始变慢；在 4×65536 时，每个条目也大约有 2.2 个周期。我相信这（笔记本电脑！）CPU 此时缓存带宽有限。
有没有办法指示 GCC（我使用的是 4.8.4）完全展开底部函数中的 while 循环 ，即剥离这个环形？循环的迭代次数在编译时已知：58.
让我先解释一下我尝试过的方法。
通过检查 GAS 输出：
gcc -fpic -O2 -S GEPDOT.c
使用了 12 个寄存器 XMM0 - XMM11。如果我将标志 -funroll-loops 传递给 gcc:
gcc -fpic -O2 -funroll-loops -S GEPDOT.c
循环只展开了两次。我检查了 GCC 优化选项。 GCC 表示 -funroll-loops 也将打开 -frename-registers，因此当 GCC 展开循环时，其寄存器分配的优先选择是使用“剩余”寄存器。但是 XMM12 - XMM15 只剩下 4 个，所以 GCC 最多只能展开 2 次。如果有 48 个而不是 16 个可用的 XMM 寄存器，GCC 将毫无问题地展开 while 循环 4 次。
然而我又做了一个实验。我首先手动展开 while 循环两次，获得一个函数 GEPDOT_2。那么在
gcc -fpic -O2 -S GEPDOT_2.c
和
gcc -fpic -O2 -funroll-loops -S GEPDOT_2.c
因为GEPDOT_2已经用完了所有寄存器，所以不执行展开。
GCC 确实会重命名寄存器以避免引入潜在的 虚假依赖。但我很确定我的 GEPDOT; 不会有这样的潜力。即使有，也不重要。我尝试自己展开循环，展开 4 次比展开 2 次快，比不展开快。当然我可以手动展开更多次，但是很繁琐。 GCC 可以为我做这个吗？谢谢。
// C file "GEPDOT.c"
#include <emmintrin.h>
void GEPDOT (double *A, double *B, double *C) {
__m128d A1_vec = _mm_load_pd(A); A += 2;
__m128d B_vec = _mm_load1_pd(B); B++;
__m128d C1_vec = A1_vec * B_vec;
__m128d A2_vec = _mm_load_pd(A); A += 2;
__m128d C2_vec = A2_vec * B_vec;
B_vec = _mm_load1_pd(B); B++;
__m128d C3_vec = A1_vec * B_vec;
__m128d C4_vec = A2_vec * B_vec;
B_vec = _mm_load1_pd(B); B++;
__m128d C5_vec = A1_vec * B_vec;
__m128d C6_vec = A2_vec * B_vec;
B_vec = _mm_load1_pd(B); B++;
__m128d C7_vec = A1_vec * B_vec;
A1_vec = _mm_load_pd(A); A += 2;
__m128d C8_vec = A2_vec * B_vec;
B_vec = _mm_load1_pd(B); B++;
int k = 58;
/* can compiler unroll the loop completely (i.e., peel this loop)? */
while (k--) {
C1_vec += A1_vec * B_vec;
A2_vec = _mm_load_pd(A); A += 2;
C2_vec += A2_vec * B_vec;
B_vec = _mm_load1_pd(B); B++;
C3_vec += A1_vec * B_vec;
C4_vec += A2_vec * B_vec;
B_vec = _mm_load1_pd(B); B++;
C5_vec += A1_vec * B_vec;
C6_vec += A2_vec * B_vec;
B_vec = _mm_load1_pd(B); B++;
C7_vec += A1_vec * B_vec;
A1_vec = _mm_load_pd(A); A += 2;
C8_vec += A2_vec * B_vec;
B_vec = _mm_load1_pd(B); B++;
}
C1_vec += A1_vec * B_vec;
A2_vec = _mm_load_pd(A);
C2_vec += A2_vec * B_vec;
B_vec = _mm_load1_pd(B); B++;
C3_vec += A1_vec * B_vec;
C4_vec += A2_vec * B_vec;
B_vec = _mm_load1_pd(B); B++;
C5_vec += A1_vec * B_vec;
C6_vec += A2_vec * B_vec;
B_vec = _mm_load1_pd(B);
C7_vec += A1_vec * B_vec;
C8_vec += A2_vec * B_vec;
/* [write-back] */
A1_vec = _mm_load_pd(C); C1_vec = A1_vec - C1_vec;
A2_vec = _mm_load_pd(C + 2); C2_vec = A2_vec - C2_vec;
A1_vec = _mm_load_pd(C + 4); C3_vec = A1_vec - C3_vec;
A2_vec = _mm_load_pd(C + 6); C4_vec = A2_vec - C4_vec;
A1_vec = _mm_load_pd(C + 8); C5_vec = A1_vec - C5_vec;
A2_vec = _mm_load_pd(C + 10); C6_vec = A2_vec - C6_vec;
A1_vec = _mm_load_pd(C + 12); C7_vec = A1_vec - C7_vec;
A2_vec = _mm_load_pd(C + 14); C8_vec = A2_vec - C8_vec;
_mm_store_pd(C,C1_vec); _mm_store_pd(C + 2,C2_vec);
_mm_store_pd(C + 4,C3_vec); _mm_store_pd(C + 6,C4_vec);
_mm_store_pd(C + 8,C5_vec); _mm_store_pd(C + 10,C6_vec);
_mm_store_pd(C + 12,C7_vec); _mm_store_pd(C + 14,C8_vec);
}
更新 1
感谢@user3386109 的评论，我想稍微扩展一下这个问题。 @user3386109 提出了一个 非常好的 问题。实际上，当有这么多并行指令要调度时，我确实怀疑编译器优化寄存器分配的能力。
我个人认为一个可靠的方法是先在 asm 内联汇编中编写循环体（这是 HPC 的关键），然后根据需要复制多次.今年早些时候我有一个不受欢迎的post：
我曾经（现在仍然）是汇编编码的婴儿，所以这对我来说是一个很好的案例研究，可以让我学习一些 x86 汇编。但是在很长的运行中，我不倾向于编码GEPDOT 占很大比例的汇编。主要有以下三个原因：
- asm 内联汇编因为不可移植而受到批评。虽然我不明白为什么。也许是因为不同的机器有不同的寄存器被破坏？
- 编译器也越来越好。所以我还是更喜欢算法优化和更好的C编码习惯来帮助编译器产生好的输出；
- 最后一个原因更重要。迭代次数可能并不总是 58。我正在开发一个高性能矩阵分解子程序。对于缓存阻塞因子
nb，迭代次数将为nb-2。我不会将nb作为函数参数，就像我在之前的 post 中所做的那样。这是一个机器特定参数，将被定义为一个宏。所以迭代次数在编译时是已知的，但可能因机器而异。猜猜我在为各种nb手动循环展开时需要做多少繁琐的工作。因此，如果有一种方法可以简单地指示编译器剥离一个循环，那就太好了。
如果您也能分享一些制作高性能但可移植的库的经验，我将不胜感激。
尝试调整优化器参数：
gcc -O3 -funroll-loops --param max-completely-peeled-insns=1000 --param max-completely-peel-times=100
这应该可以解决问题。
这不是答案，但其他尝试使用 GCC 向量化矩阵乘法的人可能会感兴趣。
下面，我假设 c 是行主顺序的 4×4 矩阵，a 是 4 行，n-列主矩阵（t运行sposed），b是一个4列， n-行矩阵行主序，运算为c = a × b + c，其中×表示矩阵乘法。
实现这个的最简单的函数是
void slow_4(double *c,
const double *a,
const double *b,
size_t n)
{
size_t row, col, i;
for (row = 0; row < 4; row++)
for (col = 0; col < 4; col++)
for (i = 0; i < n; i++)
c[4*row+col] += a[4*i+row] * b[4*i+col];
}
GCC 使用
为 SSE2/SSE3 生成了非常好的代码#if defined(__SSE2__) || defined(__SSE3__)
typedef double vec2d __attribute__((vector_size (2 * sizeof (double))));
void fast_4(vec2d *c,
const vec2d *a,
const vec2d *b,
size_t n)
{
const vec2d *const b_end = b + 2L * n;
vec2d s00 = c[0];
vec2d s02 = c[1];
vec2d s10 = c[2];
vec2d s12 = c[3];
vec2d s20 = c[4];
vec2d s22 = c[5];
vec2d s30 = c[6];
vec2d s32 = c[7];
while (b < b_end) {
const vec2d b0 = b[0];
const vec2d b2 = b[1];
const vec2d a0 = { a[0][0], a[0][0] };
const vec2d a1 = { a[0][1], a[0][1] };
const vec2d a2 = { a[1][0], a[1][0] };
const vec2d a3 = { a[1][1], a[1][1] };
s00 += a0 * b0;
s10 += a1 * b0;
s20 += a2 * b0;
s30 += a3 * b0;
s02 += a0 * b2;
s12 += a1 * b2;
s22 += a2 * b2;
s32 += a3 * b2;
b += 2;
a += 2;
}
c[0] = s00;
c[1] = s02;
c[2] = s10;
c[3] = s12;
c[4] = s20;
c[5] = s22;
c[6] = s30;
c[7] = s32;
}
#endif
对于 AVX，GCC 可以做得更好
#if defined(__AVX__) || defined(__AVX2__)
typedef double vec4d __attribute__((vector_size (4 * sizeof (double))));
void fast_4(vec4d *c,
const vec4d *a,
const vec4d *b,
size_t n)
{
const vec4d *const b_end = b + n;
vec4d s0 = c[0];
vec4d s1 = c[1];
vec4d s2 = c[2];
vec4d s3 = c[3];
while (b < b_end) {
const vec4d bc = *(b++);
const vec4d ac = *(a++);
const vec4d a0 = { ac[0], ac[0], ac[0], ac[0] };
const vec4d a1 = { ac[1], ac[1], ac[1], ac[1] };
const vec4d a2 = { ac[2], ac[2], ac[2], ac[2] };
const vec4d a3 = { ac[3], ac[3], ac[3], ac[3] };
s0 += a0 * bc;
s1 += a1 * bc;
s2 += a2 * bc;
s3 += a3 * bc;
}
c[0] = s0;
c[1] = s1;
c[2] = s2;
c[3] = s3;
}
#endif
使用gcc-4.8.4(-O2 -march=x86-64 -mtune=generic -msse3)生成的程序集的SSE3版本本质上是
fast_4:
salq , %rcx
movapd (%rdi), %xmm13
addq %rdx, %rcx
cmpq %rcx, %rdx
movapd 16(%rdi), %xmm12
movapd 32(%rdi), %xmm11
movapd 48(%rdi), %xmm10
movapd 64(%rdi), %xmm9
movapd 80(%rdi), %xmm8
movapd 96(%rdi), %xmm7
movapd 112(%rdi), %xmm6
jnb .L2
.L3:
movddup (%rsi), %xmm5
addq , %rdx
movapd -32(%rdx), %xmm1
addq , %rsi
movddup -24(%rsi), %xmm4
movapd %xmm5, %xmm14
movddup -16(%rsi), %xmm3
movddup -8(%rsi), %xmm2
mulpd %xmm1, %xmm14
movapd -16(%rdx), %xmm0
cmpq %rdx, %rcx
mulpd %xmm0, %xmm5
addpd %xmm14, %xmm13
movapd %xmm4, %xmm14
mulpd %xmm0, %xmm4
addpd %xmm5, %xmm12
mulpd %xmm1, %xmm14
addpd %xmm4, %xmm10
addpd %xmm14, %xmm11
movapd %xmm3, %xmm14
mulpd %xmm0, %xmm3
mulpd %xmm1, %xmm14
mulpd %xmm2, %xmm0
addpd %xmm3, %xmm8
mulpd %xmm2, %xmm1
addpd %xmm14, %xmm9
addpd %xmm0, %xmm6
addpd %xmm1, %xmm7
ja .L3
.L2:
movapd %xmm13, (%rdi)
movapd %xmm12, 16(%rdi)
movapd %xmm11, 32(%rdi)
movapd %xmm10, 48(%rdi)
movapd %xmm9, 64(%rdi)
movapd %xmm8, 80(%rdi)
movapd %xmm7, 96(%rdi)
movapd %xmm6, 112(%rdi)
ret
生成的程序集(-O2 -march=x86-64 -mtune=generic -mavx)的AVX版本本质上是
fast_4:
salq , %rcx
vmovapd (%rdi), %ymm5
addq %rdx, %rcx
vmovapd 32(%rdi), %ymm4
cmpq %rcx, %rdx
vmovapd 64(%rdi), %ymm3
vmovapd 96(%rdi), %ymm2
jnb .L2
.L3:
addq , %rsi
vmovapd -32(%rsi), %ymm1
addq , %rdx
vmovapd -32(%rdx), %ymm0
cmpq %rdx, %rcx
vpermilpd [=14=], %ymm1, %ymm6
vperm2f128 [=14=], %ymm6, %ymm6, %ymm6
vmulpd %ymm0, %ymm6, %ymm6
vaddpd %ymm6, %ymm5, %ymm5
vpermilpd , %ymm1, %ymm6
vperm2f128 [=14=], %ymm6, %ymm6, %ymm6
vmulpd %ymm0, %ymm6, %ymm6
vaddpd %ymm6, %ymm4, %ymm4
vpermilpd [=14=], %ymm1, %ymm6
vpermilpd , %ymm1, %ymm1
vperm2f128 , %ymm6, %ymm6, %ymm6
vperm2f128 , %ymm1, %ymm1, %ymm1
vmulpd %ymm0, %ymm6, %ymm6
vmulpd %ymm0, %ymm1, %ymm0
vaddpd %ymm6, %ymm3, %ymm3
vaddpd %ymm0, %ymm2, %ymm2
ja .L3
.L2:
vmovapd %ymm5, (%rdi)
vmovapd %ymm4, 32(%rdi)
vmovapd %ymm3, 64(%rdi)
vmovapd %ymm2, 96(%rdi)
vzeroupper
ret
寄存器调度不是最优的，我猜，但它看起来也不糟糕。我个人对上面的内容很满意，此时没有尝试手动优化它。
在 Core i5-4200U 处理器（支持 AVX2）上，上述函数的快速版本在 1843 CPU 个 SSE3 周期（中位数）和 1248 中计算两个 4×256 矩阵的乘积AVX2 的周期。这归结为每个矩阵条目有 1.8 和 1.22 个周期。未向量化的慢速版本每个矩阵条目大约需要 11 个周期，以供比较。
（循环计数是中值，即一半的测试速度更快。我只是 运行 一些粗略的基准测试，大约 100k 次重复，所以请对这些数字持保留态度。）
在这个 CPU 上，缓存效果使得 AVX2 在 4×512 矩阵大小下仍为每个条目 1.2 个周期，但在 4×1024 时，它下降到 1.4，在 4×4096 时下降到 1.5 ，在 4×8192 到 1.8，以及在 4×65536 到 2.2 每个条目的周期。 SSE3 版本保持在每个条目 1.8 个周期，最高为 4×3072，此时它开始变慢；在 4×65536 时，每个条目也大约有 2.2 个周期。我相信这（笔记本电脑！）CPU 此时缓存带宽有限。