为什么矢量化循环没有性能提升答案

【问题标题】：Why vectorizing the loop does not have performance improvement为什么矢量化循环没有性能提升
【发布时间】：2013-08-12 03:20:28
【问题描述】：

我正在研究矢量化对程序性能的影响。对此，我写了如下代码：

#include <stdio.h>
#include <sys/time.h>
#include <stdlib.h>

#define LEN 10000000

int main(){

    struct timeval stTime, endTime;

    double* a = (double*)malloc(LEN*sizeof(*a));
    double* b = (double*)malloc(LEN*sizeof(*b));
    double* c = (double*)malloc(LEN*sizeof(*c));

    int k;
    for(k = 0; k < LEN; k++){
        a[k] = rand();
        b[k] = rand();
    }

    gettimeofday(&stTime, NULL);

    for(k = 0; k < LEN; k++)
        c[k] = a[k] * b[k];

    gettimeofday(&endTime, NULL);

    FILE* fh = fopen("dump", "w");
    for(k = 0; k < LEN; k++)
        fprintf(fh, "c[%d] = %f\t", k, c[k]);
    fclose(fh);

    double timeE = (double)(endTime.tv_usec + endTime.tv_sec*1000000 - stTime.tv_usec - stTime.tv_sec*1000000);

    printf("Time elapsed: %f\n", timeE);

    return 0;
}

在这段代码中，我只是对两个向量进行初始化和相乘。结果保存在矢量c 中。我主要感兴趣的是矢量化以下循环的效果：

for(k = 0; k < LEN; k++)
    c[k] = a[k] * b[k];

我使用以下两个命令编译代码：

1) icc -O2 TestSMID.c -o TestSMID -no-vec -no-simd
2) icc -O2 TestSMID.c -o TestSMID -vec-report2

我希望看到性能改进，因为第二个命令成功地矢量化了循环。然而，我的研究表明，当循环被矢量化时，性能并没有提升。

我可能在这里遗漏了一些东西，因为我对这个话题不是很熟悉。所以，如果我的代码有问题，请告诉我。

提前感谢您的帮助。

PS：我使用的是 Mac OSX，所以不需要对齐数据，因为所有分配的内存都是 16 字节对齐的。

编辑：我首先要感谢大家的cmets和答案。我想到了@Mysticial 提出的答案，这里还有一些应该提到的点。首先，正如@Vinska 提到的，c[k]=a[k]*b[k] 不仅仅需要一个周期。除了循环索引增量和进行比较以确保k 小于LEN 之外，还需要执行其他操作来执行操作。看一下编译器生成的汇编代码，可以看出一个简单的乘法需要不止一个周期。矢量化版本如下所示：

L_B1.9:                         # Preds L_B1.8
        movq      %r13, %rax                                    #25.5
        andq      $15, %rax                                     #25.5
        testl     %eax, %eax                                    #25.5
        je        L_B1.12       # Prob 50%                      #25.5
                                # LOE rbx r12 r13 r14 r15 eax
L_B1.10:                        # Preds L_B1.9
        testb     $7, %al                                       #25.5
        jne       L_B1.32       # Prob 10%                      #25.5
                                # LOE rbx r12 r13 r14 r15
L_B1.11:                        # Preds L_B1.10
        movsd     (%r14), %xmm0                                 #26.16
        movl      $1, %eax                                      #25.5
        mulsd     (%r15), %xmm0                                 #26.23
        movsd     %xmm0, (%r13)                                 #26.9
                                # LOE rbx r12 r13 r14 r15 eax
L_B1.12:                        # Preds L_B1.11 L_B1.9
        movl      %eax, %edx                                    #25.5
        movl      %eax, %eax                                    #26.23
        negl      %edx                                          #25.5
        andl      $1, %edx                                      #25.5
        negl      %edx                                          #25.5
        addl      $10000000, %edx                               #25.5
        lea       (%r15,%rax,8), %rcx                           #26.23
        testq     $15, %rcx                                     #25.5
        je        L_B1.16       # Prob 60%                      #25.5
                                # LOE rdx rbx r12 r13 r14 r15 eax
L_B1.13:                        # Preds L_B1.12
        movl      %eax, %eax                                    #25.5
                                # LOE rax rdx rbx r12 r13 r14 r15
L_B1.14:                        # Preds L_B1.14 L_B1.13
        movups    (%r15,%rax,8), %xmm0                          #26.23
        movsd     (%r14,%rax,8), %xmm1                          #26.16
        movhpd    8(%r14,%rax,8), %xmm1                         #26.16
        mulpd     %xmm0, %xmm1                                  #26.23
        movntpd   %xmm1, (%r13,%rax,8)                          #26.9
        addq      $2, %rax                                      #25.5
        cmpq      %rdx, %rax                                    #25.5
        jb        L_B1.14       # Prob 99%                      #25.5
        jmp       L_B1.20       # Prob 100%                     #25.5
                                # LOE rax rdx rbx r12 r13 r14 r15
L_B1.16:                        # Preds L_B1.12
        movl      %eax, %eax                                    #25.5
                                # LOE rax rdx rbx r12 r13 r14 r15
L_B1.17:                        # Preds L_B1.17 L_B1.16
        movsd     (%r14,%rax,8), %xmm0                          #26.16
        movhpd    8(%r14,%rax,8), %xmm0                         #26.16
        mulpd     (%r15,%rax,8), %xmm0                          #26.23
        movntpd   %xmm0, (%r13,%rax,8)                          #26.9
        addq      $2, %rax                                      #25.5
        cmpq      %rdx, %rax                                    #25.5
        jb        L_B1.17       # Prob 99%                      #25.5
                                # LOE rax rdx rbx r12 r13 r14 r15
L_B1.18:                        # Preds L_B1.17
        mfence                                                  #25.5
                                # LOE rdx rbx r12 r13 r14 r15
L_B1.19:                        # Preds L_B1.18
        mfence                                                  #25.5
                                # LOE rdx rbx r12 r13 r14 r15
L_B1.20:                        # Preds L_B1.14 L_B1.19 L_B1.32
        cmpq      $10000000, %rdx                               #25.5
        jae       L_B1.24       # Prob 0%                       #25.5
                                # LOE rdx rbx r12 r13 r14 r15
L_B1.22:                        # Preds L_B1.20 L_B1.22
        movsd     (%r14,%rdx,8), %xmm0                          #26.16
        mulsd     (%r15,%rdx,8), %xmm0                          #26.23
        movsd     %xmm0, (%r13,%rdx,8)                          #26.9
        incq      %rdx                                          #25.5
        cmpq      $10000000, %rdx                               #25.5
        jb        L_B1.22       # Prob 99%                      #25.5
                                # LOE rdx rbx r12 r13 r14 r15
L_B1.24:                        # Preds L_B1.22 L_B1.20

非矢量化版本是：

L_B1.9:                         # Preds L_B1.8
        xorl      %eax, %eax                                    #25.5
                                # LOE rbx r12 r13 r14 r15 eax
L_B1.10:                        # Preds L_B1.10 L_B1.9
        lea       (%rax,%rax), %edx                             #26.9
        incl      %eax                                          #25.5
        cmpl      $5000000, %eax                                #25.5
        movsd     (%r15,%rdx,8), %xmm0                          #26.16
        movsd     8(%r15,%rdx,8), %xmm1                         #26.16
        mulsd     (%r13,%rdx,8), %xmm0                          #26.23
        mulsd     8(%r13,%rdx,8), %xmm1                         #26.23
        movsd     %xmm0, (%rbx,%rdx,8)                          #26.9
        movsd     %xmm1, 8(%rbx,%rdx,8)                         #26.9
        jb        L_B1.10       # Prob 99%                      #25.5
                                # LOE rbx r12 r13 r14 r15 eax

除此之外，处理器不仅仅加载 24 个字节。在每次访问内存时，都会加载一个完整的行（64 字节）。更重要的是，由于a、b 和c 所需的内存是连续的，预取器肯定会有很大帮助并提前加载下一个块。说了这么多，我觉得@Mysticial计算的内存带宽太悲观了。

此外，Intel Vectorization Guide 中提到了使用 SIMD 来提高程序性能的非常简单的加法。因此，对于这个非常简单的循环，我们似乎应该能够获得一些性能提升。

编辑2：再次感谢您的 cmets。另外，感谢@Mysticial 示例代码，我终于看到了SIMD 对性能提升的影响。正如 Mysticial 所提到的，问题在于内存带宽。通过为适合L1缓存的a、b和c选择小尺寸，可以看出SIMD有助于显着提高性能。以下是我得到的结果：

icc -O2 -o TestSMIDNoVec -no-vec TestSMID2.c: 17.34 sec

icc -O2 -o TestSMIDVecNoUnroll -vec-report2 TestSMID2.c: 9.33 sec

展开循环可以进一步提高性能：

icc -O2 -o TestSMIDVecUnroll -vec-report2 TestSMID2.c -unroll=8: 8.6sec

另外，我应该提一下，当使用-O2 编译时，我的处理器只需要一个周期即可完成一次迭代。

PS：我的电脑是 Macbook Pro core i5 @2.5GHz（双核）

【问题讨论】：

我刚刚更新了我的答案，以证明我的处理器每个周期能够进行 1 次迭代，并解释了它是如何实现的。
我真的很讨厌提出这个问题，但是构建命令会将两个版本的可执行文件放在同一个文件中。如果两个版本有不同的名称会更清楚。
您说“不需要对齐”，但是生成的 asm 代码会检查所有对齐的可能性。有一个未对齐的 srces 循环，以及一个使用带有内存操作数的 mulpd 的循环。然而，即使是对齐版本也使用奇怪的movsd + movhpd 序列来加载 128b。我认为这适用于 c 和 a 对齐，b 未对齐（在标量介绍之后）。我想我记得在一些较旧的架构上读到过，2 insn 序列有时比movupd 快。循环的唯一目标对齐版本对一个源使用 movupd，对另一个源使用 2 insn 方法，/boggle。
LEN你选了什么尺寸的？

标签： c performance simd icc

【解决方案1】：

这个原始答案早在 2013 年就有效。截至 2017 年硬件，情况发生了很大变化，问题和答案都已过时。

2017 年更新请参阅此答案的结尾。

原答案（2013 年）：

因为您受到内存带宽的限制。

虽然矢量化和其他微优化可以提高计算速度，但它们无法提高内存速度。

在你的例子中：

for(k = 0; k < LEN; k++)
    c[k] = a[k] * b[k];

你只做很少的工作就对所有内存进行了一次遍历。这会耗尽您的内存带宽。

因此，无论如何优化（矢量化、展开等），它都不会变得更快。

2013 年的典型台式机具有大约 10 GB/s 的内存带宽*。
您的循环涉及 24 字节/迭代。 p>

如果没有矢量化，现代 x64 处理器可能每个周期可以进行大约 1 次迭代*。

假设您以 4 GHz 运行：

(4 * 10^9) * 24 bytes/iteration = 96 GB/s

这几乎是内存带宽的 10 倍 - 没有矢量化。

*毫不奇怪，一些人怀疑我上面给出的数字，因为我没有给出引用。好吧，从经验来看，这些都是我的想法。所以这里有一些基准来证明这一点。

循环迭代可以以 1 个循环/迭代的速度运行：

如果我们减少 LEN 使其适合缓存，我们可以摆脱内存瓶颈。
（我在 C++ 中对此进行了测试，因为它更容易。但这没有区别。）

#include <iostream>
#include <time.h>
using std::cout;
using std::endl;

int main(){
    const int LEN = 256;

    double *a = (double*)malloc(LEN*sizeof(*a));
    double *b = (double*)malloc(LEN*sizeof(*a));
    double *c = (double*)malloc(LEN*sizeof(*a));

    int k;
    for(k = 0; k < LEN; k++){
        a[k] = rand();
        b[k] = rand();
    }

    clock_t time0 = clock();

    for (int i = 0; i < 100000000; i++){
        for(k = 0; k < LEN; k++)
            c[k] = a[k] * b[k];
    }

    clock_t time1 = clock();
    cout << (double)(time1 - time0) / CLOCKS_PER_SEC << endl;
}

处理器：Intel Core i7 2600K @ 4.2 GHz
编译器：Visual Studio 2012
时间：6.55 秒

在此测试中，我仅在 6.55 秒内运行了 25,600,000,000 次迭代。

6.55 * 4.2 GHz = 27,510,000,000 个周期
27,510,000,000 / 25,600,000,000 = 1.074 个周期/迭代

现在，如果您想知道如何做到这一点：

2 次加载
1 家商店
1 倍
增量计数器
比较+分支

一个循环...

这是因为现代处理器和编译器非常棒。

虽然每个操作都有延迟（尤其是乘法），但处理器能够同时执行多个迭代。我的测试机器是 Sandy Bridge 处理器，它能够承受 2x128b 负载、1x128b 存储和 1x256b 矢量 FP 乘以每个周期。如果负载是微融合微指令的内存源操作数，则可能还有另外一两个向量或整数运算。（2 次加载 + 1 次存储吞吐量仅在使用 256b AVX 加载/存储时，否则每个周期只有两个总内存操作（最多一个存储）。

查看程序集（为简洁起见，我将省略），编译器似乎展开了循环，从而减少了循环开销。但它并没有很好地对其进行矢量化。

内存带宽约为 10 GB/s：

最简单的测试方法是通过memset()：

#include <iostream>
#include <time.h>
using std::cout;
using std::endl;

int main(){
    const int LEN = 1 << 30;    //  1GB

    char *a = (char*)calloc(LEN,1);

    clock_t time0 = clock();

    for (int i = 0; i < 100; i++){
        memset(a,0xff,LEN);
    }

    clock_t time1 = clock();
    cout << (double)(time1 - time0) / CLOCKS_PER_SEC << endl;
}

处理器：Intel Core i7 2600K @ 4.2 GHz
编译器：Visual Studio 2012
时间：5.811 秒

所以我的机器需要 5.811 秒才能写入 100 GB 内存。这大约是 17.2 GB/s。

而且我的处理器更高端。 Nehalem 和 Core 2 代处理器的内存带宽更少。

2017 年 3 月更新：

截至 2017 年，情况变得更加复杂。

多亏了 DDR4 和四通道内存，单线程不再可能使内存带宽饱和。但带宽问题并不一定会消失。尽管带宽增加了，但处理器内核也得到了改进——而且数量更多。

从数学上来说：

每个内核都有带宽限制X。
主内存的带宽限制为Y。
在旧系统上，X > Y。
在当前的高端系统上，X < Y。但是X * (# of cores) > Y。

早在 2013 年：Sandy Bridge @ 4 GHz + 双通道 DDR3 @ 1333 MHz

无矢量化（8 字节加载/存储）：X = 32 GB/s 和 Y = ~17 GB/s
矢量化 SSE*（16 字节加载/存储）：X = 64 GB/s 和 Y = ~17 GB/s

现在在 2017 年：Haswell-E @ 4 GHz + 四通道 DDR4 @ 2400 MHz

无矢量化（8 字节加载/存储）：X = 32 GB/s 和 Y = ~70 GB/s
矢量化 AVX*（32 字节加载/存储）：X = 64 GB/s 和 Y = ~70 GB/s

_{（对于 Sandy Bridge 和 Haswell，无论 SIMD 宽度如何，缓存中的架构限制都会将带宽限制为大约 16 字节/周期。）}

所以现在，单个线程并不总是能够使内存带宽饱和。您需要进行矢量化以达到X 的限制。但是您仍然会在 2 个或更多线程时达到 Y 的主内存带宽限制。

但有一件事并没有改变，而且很可能在很长一段时间内都不会改变：如果总内存带宽不饱和，您将无法在所有内核上运行占用带宽的循环。

【讨论】：

感谢您的回答。你说的对。我把事情复杂化了，体验了性能的提升。
+1：这需要出现在常见问题解答中或成为“首选”答案 - 大部分初学者优化问题似乎都属于这一类。
如果我们用-O0编译它会怎样？ CPU 是否在一个周期内执行每次迭代？
@matmul 它仅在您重用数据时有效。如果所有东西都只被触摸一次，那么就没有什么可以做的了。
@Zboson 显然这取决于机器。在具有多个 NUMA 节点的机器上，您不太可能在单线程上获得全部带宽。在 Haswell-E 上，内存足够快，您可能需要矢量化以仅使用单个线程来最大化带宽。话虽如此，但这并没有带走重点。这个问题的代码迟早会遇到带宽问题。

【解决方案2】：

正如 Mysticial 已经描述的，主内存带宽限制是大缓冲区的瓶颈。解决这个问题的方法是重新设计您的处理以在适合缓存的块中工作。（而不是将整个 200MiB 的双精度相乘，而是仅相乘 128kiB，然后对其进行处理。因此，使用乘法输出的代码会发现它仍在 L2 缓存中。L2 通常为 256kiB，并且是每个 CPU 内核私有的，关于最近的英特尔设计。）

这种技术称为cache blocking 或loop tiling。对于某些算法来说可能会很棘手，但收益是二级缓存带宽与主内存带宽之间的差异。

如果您这样做，请确保编译器未仍在生成流式存储 (movnt...)。这些写入绕过缓存以避免不适合的数据污染它。下一次读取该数据需要触及主存。

【讨论】：

【解决方案3】：

编辑：修改了答案很多。另外，请忽略我之前写的关于 Mystical 的答案并不完全正确的大部分内容。尽管如此，我仍然不同意它受到内存的限制，因为尽管进行了各种各样的测试，但我看不到任何原始代码受内存速度限制的迹象。与此同时，它不断显示出受 CPU 限制的明显迹象。

可能有很多原因。而且由于原因可能非常依赖于硬件，因此我决定不应该根据猜测进行推测。只是概述一下我在以后的测试中遇到的这些事情，我使用了一种更准确可靠的 CPU 时间测量方法和循环 1000 次。我相信这些信息可能会有所帮助。但请谨慎对待，因为它取决于硬件。

使用 SSE 系列指令时，我得到的矢量化代码比非矢量化代码快 10% 以上。
使用 SSE 系列的矢量化代码和使用 AVX 的矢量化代码运行时或多或少具有相同的性能。
使用 AVX 指令时，非矢量化代码运行速度最快 - 比我尝试过的所有其他方法快 25% 或更多。
在所有情况下，结果都与 CPU 时钟成线性关系。
结果几乎不受内存时钟的影响。
内存延迟对结果的影响很大 - 远大于内存时钟，但不如 CPU 时钟对结果的影响那么大。

WRT Mystical 的每个时钟运行近 1 次迭代的示例 - 我没想到 CPU 调度程序会如此高效，并且假设每 1.5-2 个时钟周期进行 1 次迭代。但令我惊讶的是，事实并非如此。我确实错了，对此感到抱歉。我自己的 CPU 运行它的效率更高 - 1.048 个周期/迭代。所以我可以证明 Mystical 的这部分答案是绝对正确的。

【讨论】：

Along with the multiply instruction, the code of the loop has to execute several other instructions as well, including the conditional 啊，你没有向我们展示真实代码。在循环中添加条件将有效地破坏分支预测。顺便说一句，您报告的百分之几的收益是徒劳的。您仍然受到总线带宽的约束。恕我直言，手动展开只会导致更少的分支预测未命中，因为迭代次数更少。 L1局部性基本相同。
@wildplasser 定义“真实代码”。还有其他一些事情：数据的总大小为 10,000,000 * 8 * 3 = 228 MB。在我的正常时钟上，我的理论内存带宽是 29.8 GB/s。如果我将 CPU 设置为可用的最低时钟速度，这部分代码将运行大约 1.1 秒。在那段时间内，它可以发送整个数据 131 次。所以我看不到内存瓶颈会发生在哪里。此外，“内存瓶颈”理论与以下事实不符，即如果我将 CPU 时钟加倍，那部分代码的运行速度就会提高一倍，而将内存时钟加倍几乎没有任何作用。
@wildplasser 另外，百分之几？最快的非矢量化和最快的矢量化之间的差异略高于 6.5%。这可能看起来不多，但在更大的范围内可能非常重要。有了这样的差异，这将意味着例如花费 11 小时 20 分钟的 CPU 时间，而不是花费 12 小时。 40 分钟。小事情确实加起来，所以它远非“徒劳”
复制到自动存储可以避免/减少 L2 缓存的影响，这里减少了 30%。我将其添加为答案，因为我需要格式化。
WRT real code ：我首先以为你是 OP。对不起！

【解决方案4】：

以防万一 a[] b[] 和 c[] 争夺 L2 缓存 ::

#include <string.h> /* for memcpy */

 ...

 gettimeofday(&stTime, NULL);

    for(k = 0; k < LEN; k += 4) {
        double a4[4], b4[4], c4[4];
        memcpy(a4,a+k, sizeof a4);
        memcpy(b4,b+k, sizeof b4);
        c4[0] = a4[0] * b4[0];
        c4[1] = a4[1] * b4[1];
        c4[2] = a4[2] * b4[2];
        c4[3] = a4[3] * b4[3];
        memcpy(c+k,c4, sizeof c4);
        }

    gettimeofday(&endTime, NULL);

将运行时间从 98429.000000 减少到 67213.000000；展开循环 8 倍将其减少到 57157.000000。

【讨论】：

对我来说，它的值要小得多——仅比 OP 的香草版本增加了 2%。（4 折和 8 折展开的结果相同）
当我进行优化时，我的收获消失了。 GCC 似乎会自动展开循环，而且它似乎也在以某种方式按摩缓存。