为什么将循环向量化到64位元素上，对于大缓冲区没有性能提升？

Question

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

#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字节对齐的。

编辑： 首先，我想感谢大家的评论和答案。 我考虑了@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个字节）。更重要的是，由于所需的内存是连续的，预取器肯定会有很大帮助，并提前加载下一个块。

话虽如此，我认为@Mysticial计算的内存带宽过于悲观。

此外，在Intel向量化指南中提到，使用SIMD来改善非常简单的加法程序的性能。因此，似乎我们应该能够在这个非常简单的循环中获得一些性能提升。

编辑2： 再次感谢您的评论。还要感谢@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核心i5 @ 2.5GHz（双核）

- Pouya

---

我刚刚更新了我的答案，证明我的处理器能够每个周期执行1次迭代，并解释了这是如何可能的。 - Mysticial

我真的不想提这个问题，但是构建命令会将两个可执行文件版本放在同一个文件中。如果这两个版本有不同的名称，那么会更加清晰明了。 - wallyk

你说“不需要对齐”，但生成的汇编代码检查了所有对齐可能性。有一个循环用于未对齐的源，另一个使用带有内存操作数的mulpd。然而，即使是对齐版本也使用奇怪的movsd+movhpd序列来加载128b。我认为这是针对c和a对齐，b未对齐（在标量介绍后）的情况。我记得在一些旧架构上，2个指令序列有时比movupd更快。只有目标对齐版本的循环使用movupd作为一个源，另一个源使用2个指令方法，/boggle。 - Peter Cordes

你选择了什么样的 LEN 大小？ - Manolete

Answer 1

这个原始答案在2013年是正确的。但随着2017年硬件的变化，问题和答案都已经过时了。

请参考本答案末尾的2017年更新。

---

原始答案（2013）：

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

虽然矢量化和其他微小优化可以提高计算速度，但它们不能增加内存速度。

以您的示例为例：

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

您在所有内存上进行了一次单次遍历，做的工作非常少。这使得您的内存带宽达到了最大值。

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

2013年的典型台式机具有约10 GB/s的内存带宽。您的循环每次迭代使用24字节。

没有矢量化的情况下，现代x64处理器可能每个周期执行1次迭代。

假设您运行速度为4 GHz：

(4 * 10 ^ 9) * 24字节/迭代=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秒就完成了256亿次迭代。

• 6.55 * 4.2 GHz = 27,510,000,000 cycles
• 27,510,000,000 / 25,600,000,000 = 1.074 cycles/iteration

---

如果你想知道如何做到以下操作：

• 2次加载
• 1次存储
• 1次乘法
• 增加计数器
• 比较+跳转

所有这些操作都在一个时钟周期内完成……

那是因为现代处理器和编译器非常强大。

虽然每个操作都有延迟（特别是乘法），但是处理器能够同时执行多个迭代。我的测试机器是Sandy Bridge处理器，可以每个时钟周期内持续执行2x128b的加载、1x128b的存储和1x256b的向量浮点乘法。如果加载是微操作的内存源操作数，则还可以进行另外一到两个向量或整数操作。 (使用256b AVX加载/存储时，2个加载+1个存储吞吐量，否则每个周期最多只有两个内存操作（至多一个存储）)。

从汇编代码（为简洁起见省略）来看，编译器展开了循环，从而减少了循环开销。但它并没有完全将其向量化。

---

内存带宽大约为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 * (# of cores) > Y。但是X < Y。

回到2013年：Sandy Bridge @ 4 GHz + dual-channel 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 + quad-channel 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，缓存中的架构限制将限制带宽约为每个周期16字节，而不管SIMD宽度如何。）}

因此，现在单个线程并不总是能够饱和内存带宽。您需要进行向量化才能实现X的限制。但是，使用2个或更多线程仍将达到主内存带宽限制Y。

但是有一件事情没有改变，而且可能在很长一段时间内都不会改变：如果您在所有核心上运行带宽占用率很高的循环，那么总内存带宽就会饱和。

- Mysticial

---

谢谢你的回答。你是对的。我把事情复杂化了，但也体验到了性能的提升。 - Pouya

这段话需要翻译为：+1：这应该成为常见问题解答或首选答案，因为很多初学者的优化问题似乎都属于这个类别。 - Paul R

只有在重复使用数据时，@matmul 才能起作用。如果每个数据值只被访问一次，那么就不会有太多可以做的了。 - Mysticial

你的说法“速度不会再快了”并不一定准确。特别是在你测试的机器上，我认为使用多线程可以获得至少50%以上的带宽。 - Z boson

显然这取决于设备。在具有多个NUMA节点的机器上，单线程不太可能获得完整的带宽。在Haswell-E上，内存速度足够快，您可能需要矢量化才能通过单个线程达到最大带宽。尽管如此，这并不影响要点。该问题中的代码迟早会遇到带宽问题。 - Mysticial

Answer 2

正如Mysticial所描述的那样，主存带宽限制是大缓冲区的瓶颈。解决这个问题的方法是重新设计处理方式，使其适合缓存中的块。(不要乘以整个200MiB的双精度浮点数，只需乘以128kiB，然后对其进行处理。这样使用乘积输出的代码将在L2高速缓存中找到它。L2通常为256kiB，并且在最近的英特尔设计中是每个CPU核心私有的。) 这种技术称为缓存阻塞或循环分块。 对于某些算法可能会有些棘手，但回报是L2缓存带宽与主存带宽之间的差异。

如果您这样做，请确保编译器没有仍在生成流式存储(movnt...)。这些写操作绕过缓存，以避免污染无法容纳的数据。下一个读取该数据的操作将需要访问主内存。

- Peter Cordes

Answer 3

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

---

可能有很多原因。由于原因可能与硬件相关，我决定不根据猜测进行推测。 只是要概述我在后期测试中遇到的这些事情，其中我使用了更准确可靠的CPU时间测量方法和循环1000次。我相信这些信息可能会有所帮助。但请将其视为硬件相关的一种参考。

• 使用SSE系列指令时，我得到的矢量化代码比非矢量化代码快10%以上。
• 使用SSE系列的矢量化代码和使用AVX的矢量化代码的性能基本相同。
• 使用AVX指令时，非矢量化代码运行最快-比我尝试的其他任何东西都快25%或更多。
• 结果在所有情况下都与CPU时钟线性缩放。
• 结果几乎不受内存时钟的影响。
• 结果受内存延迟的影响很大-比内存时钟影响的要多得多，但远远不及CPU时钟影响结果的程度。

关于Mystical运行近1个时钟周期的示例-我没有预料到CPU调度程序会那么有效率，并且假设每1.5-2个时钟节拍运行1次。但令我惊讶的是，事实并非如此；我肯定是错了，对此感到抱歉。我的CPU甚至更有效地运行它-每次迭代1.048个周期。因此，我可以证明Mystical答案中的这部分绝对正确。

- librin.so.1

---

除了乘法指令，循环的代码还必须执行其他几个指令，包括条件语句。啊，你没有给我们展示真正的代码。在循环内添加条件语句将有效地破坏分支预测。顺便说一下，你报告的少数收益是徒劳的。你仍然受到总线带宽的限制。在我看来，手动展开只会导致更少的分支预测失误，因为迭代次数更少。L1局部性基本相同。 - wildplasser

@wildplasser 定义“真正的代码”。还有一些其他的事情：数据的总大小为10,000,000 * 8 * 3 = 228兆字节。在我的普通时钟下，我的理论内存带宽为29.8 GB/s。如果我将CPU设置为最低可用时钟速度，则该代码部分运行约1.1秒。在那段时间里，它可以发送整个数据131次。因此，我不知道哪里会出现内存瓶颈。此外，“内存瓶颈”理论与事实不符，即如果我将CPU时钟加倍，该代码部分的运行速度将加倍，而将内存时钟加倍几乎没有任何作用。 - librin.so.1

@wildplasser，还有，几个百分点？最快的非矢量化和最快的矢量化之间的差异略大于6.5％。这可能看起来不是很多，但在更大的规模上可能非常重要。有了这样的差异，这意味着例如花费11小时20分钟的CPU时间而不是花费12小时。高达40分钟。小事情确实会累积，所以远非“徒劳无功”。 - librin.so.1

将数据复制到自动存储器可以避免/减少L2缓存效应，这里可以节省30％。我会将其添加为答案，因为我需要格式。 - wildplasser

关于“真实代码”：我最初以为你是楼主。抱歉！ - wildplasser

顺便提一下，我更新了帖子的第一部分。除其他内容外，我还包括了指令计数。 - librin.so.1

Answer 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。

- wildplasser

---

对我来说，这只是比 OP 的基本版本略高 2% 的微小增长。（使用 4 倍和 8 倍展开都得到相同的结果） - librin.so.1

当我打开优化时，我的收益消失了。GCC似乎会自动展开循环，并以某种方式调整缓存。 - wildplasser