为什么这个 SIMD 乘法不比非 SIMD 乘法快？

Question

假设我们有一个函数将两个数组相乘，每个数组有 1000000 个双精度数。在 C/C++ 中，函数如下所示：

void mul_c(double* a, double* b)
{
    for (int i = 0; i != 1000000; ++i)
    {
        a[i] = a[i] * b[i];
    }
}

编译器使用-O2 生成以下程序集：

mul_c(double*, double*):
        xor     eax, eax
.L2:
        movsd   xmm0, QWORD PTR [rdi+rax]
        mulsd   xmm0, QWORD PTR [rsi+rax]
        movsd   QWORD PTR [rdi+rax], xmm0
        add     rax, 8
        cmp     rax, 8000000
        jne     .L2
        rep ret

从上面的程序集中，编译器似乎使用了 SIMD 指令，但它每次迭代只乘以一倍。所以我决定改为在内联汇编中编写相同的函数，在那里我充分利用xmm0 寄存器并一次将两个双精度数相乘：

void mul_asm(double* a, double* b)
{
    asm volatile
    (
        ".intel_syntax noprefix             \n\t"
        "xor    rax, rax                    \n\t"
        "0:                                 \n\t"
        "movupd xmm0, xmmword ptr [rdi+rax] \n\t"
        "mulpd  xmm0, xmmword ptr [rsi+rax] \n\t"
        "movupd xmmword ptr [rdi+rax], xmm0 \n\t"
        "add    rax, 16                     \n\t"
        "cmp    rax, 8000000                \n\t"
        "jne    0b                          \n\t"
        ".att_syntax noprefix               \n\t"

        : 
        : "D" (a), "S" (b)
        : "memory", "cc"
    );
}

分别测量这两个函数的执行时间后，似乎它们都需要 1 ms 才能完成：

> gcc -O2 main.cpp
> ./a.out < input

mul_c: 1 ms
mul_asm: 1 ms

[a lot of doubles...]

我希望 SIMD 实现的速度至少是乘法/内存指令数量的一半（0 ms）的两倍。

所以我的问题是：当 SIMD 实现只执行一半的乘法/内存指令时，为什么 SIMD 实现不比普通 C/C++ 实现快？

这是完整的程序：

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

void mul_c(double* a, double* b)
{
    for (int i = 0; i != 1000000; ++i)
    {
        a[i] = a[i] * b[i];
    }
}

void mul_asm(double* a, double* b)
{
    asm volatile
    (
        ".intel_syntax noprefix             \n\t"
        "xor    rax, rax                    \n\t"
        "0:                                 \n\t"
        "movupd xmm0, xmmword ptr [rdi+rax] \n\t"
        "mulpd  xmm0, xmmword ptr [rsi+rax] \n\t"
        "movupd xmmword ptr [rdi+rax], xmm0 \n\t"
        "add    rax, 16                     \n\t"
        "cmp    rax, 8000000                \n\t"
        "jne    0b                          \n\t"
        ".att_syntax noprefix               \n\t"

        : 
        : "D" (a), "S" (b)
        : "memory", "cc"
    );
}

int main()
{
    struct timeval t1;
    struct timeval t2;
    unsigned long long time;

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

    for (int i = 0; i != 1000000; ++i)
    {
        double v;
        scanf("%lf", &v);
        a[i] = v;
        b[i] = v;
        c[i] = v;
    }

    gettimeofday(&t1, NULL);
    mul_c(a, b);
    gettimeofday(&t2, NULL);
    time = 1000 * (t2.tv_sec - t1.tv_sec) + (t2.tv_usec - t1.tv_usec) / 1000;
    printf("mul_c: %llu ms\n", time);

    gettimeofday(&t1, NULL);
    mul_asm(b, c);
    gettimeofday(&t2, NULL);
    time = 1000 * (t2.tv_sec - t1.tv_sec) + (t2.tv_usec - t1.tv_usec) / 1000;
    printf("mul_asm: %llu ms\n\n", time);

    for (int i = 0; i != 1000000; ++i)
    {
        printf("%lf\t\t\t%lf\n", a[i], b[i]);
    }

    return 0;
}

我还尝试利用所有xmm 寄存器（0-7）并删除指令依赖以获得更好的并行计算：

void mul_asm(double* a, double* b)
{
    asm volatile
    (
        ".intel_syntax noprefix                 \n\t"
        "xor    rax, rax                        \n\t"
        "0:                                     \n\t"
        "movupd xmm0, xmmword ptr [rdi+rax]     \n\t"
        "movupd xmm1, xmmword ptr [rdi+rax+16]  \n\t"
        "movupd xmm2, xmmword ptr [rdi+rax+32]  \n\t"
        "movupd xmm3, xmmword ptr [rdi+rax+48]  \n\t"
        "movupd xmm4, xmmword ptr [rdi+rax+64]  \n\t"
        "movupd xmm5, xmmword ptr [rdi+rax+80]  \n\t"
        "movupd xmm6, xmmword ptr [rdi+rax+96]  \n\t"
        "movupd xmm7, xmmword ptr [rdi+rax+112] \n\t"
        "mulpd  xmm0, xmmword ptr [rsi+rax]     \n\t"
        "mulpd  xmm1, xmmword ptr [rsi+rax+16]  \n\t"
        "mulpd  xmm2, xmmword ptr [rsi+rax+32]  \n\t"
        "mulpd  xmm3, xmmword ptr [rsi+rax+48]  \n\t"
        "mulpd  xmm4, xmmword ptr [rsi+rax+64]  \n\t"
        "mulpd  xmm5, xmmword ptr [rsi+rax+80]  \n\t"
        "mulpd  xmm6, xmmword ptr [rsi+rax+96]  \n\t"
        "mulpd  xmm7, xmmword ptr [rsi+rax+112] \n\t"
        "movupd xmmword ptr [rdi+rax], xmm0     \n\t"
        "movupd xmmword ptr [rdi+rax+16], xmm1  \n\t"
        "movupd xmmword ptr [rdi+rax+32], xmm2  \n\t"
        "movupd xmmword ptr [rdi+rax+48], xmm3  \n\t"
        "movupd xmmword ptr [rdi+rax+64], xmm4  \n\t"
        "movupd xmmword ptr [rdi+rax+80], xmm5  \n\t"
        "movupd xmmword ptr [rdi+rax+96], xmm6  \n\t"
        "movupd xmmword ptr [rdi+rax+112], xmm7 \n\t"
        "add    rax, 128                        \n\t"
        "cmp    rax, 8000000                    \n\t"
        "jne    0b                              \n\t"
        ".att_syntax noprefix                   \n\t"

        : 
        : "D" (a), "S" (b)
        : "memory", "cc"
    );
}

但它仍然以 1 毫秒的速度运行，与普通 C/C++ 实现的速度相同。

---

更新

正如 answers/cmets 所建议的，我已经实现了另一种测量执行时间的方法：

#include <stdio.h>
#include <stdlib.h>

void mul_c(double* a, double* b)
{
    for (int i = 0; i != 1000000; ++i)
    {
        a[i] = a[i] * b[i];
    }
}

void mul_asm(double* a, double* b)
{
    asm volatile
    (
        ".intel_syntax noprefix             \n\t"
        "xor    rax, rax                    \n\t"
        "0:                                 \n\t"
        "movupd xmm0, xmmword ptr [rdi+rax] \n\t"
        "mulpd  xmm0, xmmword ptr [rsi+rax] \n\t"
        "movupd xmmword ptr [rdi+rax], xmm0 \n\t"
        "add    rax, 16                     \n\t"
        "cmp    rax, 8000000                \n\t"
        "jne    0b                          \n\t"
        ".att_syntax noprefix               \n\t"

        : 
        : "D" (a), "S" (b)
        : "memory", "cc"
    );
}

void mul_asm2(double* a, double* b)
{
    asm volatile
    (
        ".intel_syntax noprefix                 \n\t"
        "xor    rax, rax                        \n\t"
        "0:                                     \n\t"
        "movupd xmm0, xmmword ptr [rdi+rax]     \n\t"
        "movupd xmm1, xmmword ptr [rdi+rax+16]  \n\t"
        "movupd xmm2, xmmword ptr [rdi+rax+32]  \n\t"
        "movupd xmm3, xmmword ptr [rdi+rax+48]  \n\t"
        "movupd xmm4, xmmword ptr [rdi+rax+64]  \n\t"
        "movupd xmm5, xmmword ptr [rdi+rax+80]  \n\t"
        "movupd xmm6, xmmword ptr [rdi+rax+96]  \n\t"
        "movupd xmm7, xmmword ptr [rdi+rax+112] \n\t"
        "mulpd  xmm0, xmmword ptr [rsi+rax]     \n\t"
        "mulpd  xmm1, xmmword ptr [rsi+rax+16]  \n\t"
        "mulpd  xmm2, xmmword ptr [rsi+rax+32]  \n\t"
        "mulpd  xmm3, xmmword ptr [rsi+rax+48]  \n\t"
        "mulpd  xmm4, xmmword ptr [rsi+rax+64]  \n\t"
        "mulpd  xmm5, xmmword ptr [rsi+rax+80]  \n\t"
        "mulpd  xmm6, xmmword ptr [rsi+rax+96]  \n\t"
        "mulpd  xmm7, xmmword ptr [rsi+rax+112] \n\t"
        "movupd xmmword ptr [rdi+rax], xmm0     \n\t"
        "movupd xmmword ptr [rdi+rax+16], xmm1  \n\t"
        "movupd xmmword ptr [rdi+rax+32], xmm2  \n\t"
        "movupd xmmword ptr [rdi+rax+48], xmm3  \n\t"
        "movupd xmmword ptr [rdi+rax+64], xmm4  \n\t"
        "movupd xmmword ptr [rdi+rax+80], xmm5  \n\t"
        "movupd xmmword ptr [rdi+rax+96], xmm6  \n\t"
        "movupd xmmword ptr [rdi+rax+112], xmm7 \n\t"
        "add    rax, 128                        \n\t"
        "cmp    rax, 8000000                    \n\t"
        "jne    0b                              \n\t"
        ".att_syntax noprefix                   \n\t"

        : 
        : "D" (a), "S" (b)
        : "memory", "cc"
    );
}

unsigned long timestamp()
{
    unsigned long a;

    asm volatile
    (
        ".intel_syntax noprefix \n\t"
        "xor   rax, rax         \n\t"
        "xor   rdx, rdx         \n\t"
        "RDTSCP                 \n\t"
        "shl   rdx, 32          \n\t"
        "or    rax, rdx         \n\t"
        ".att_syntax noprefix   \n\t"

        : "=a" (a)
        : 
        : "memory", "cc"
    );

    return a;
}

int main()
{
    unsigned long t1;
    unsigned long t2;

    double* a;
    double* b;

    a = (double*)malloc(sizeof(double) * 1000000);
    b = (double*)malloc(sizeof(double) * 1000000);

    for (int i = 0; i != 1000000; ++i)
    {
        double v;
        scanf("%lf", &v);
        a[i] = v;
        b[i] = v;
    }

    t1 = timestamp();
    mul_c(a, b);
    //mul_asm(a, b);
    //mul_asm2(a, b);
    t2 = timestamp();
    printf("mul_c: %lu cycles\n\n", t2 - t1);

    for (int i = 0; i != 1000000; ++i)
    {
        printf("%lf\t\t\t%lf\n", a[i], b[i]);
    }

    return 0;
}

当我用这个测量值运行程序时，我得到了这个结果：

mul_c:    ~2163971628 cycles
mul_asm:  ~2532045184 cycles
mul_asm2: ~5230488    cycles <-- what???

这里有两点值得注意，首先，周期数变化很大，我认为这是因为操作系统允许其他进程在其间运行。有什么方法可以防止这种情况发生，或者只在我的程序执行时计算周期？此外，mul_asm2 产生的输出与其他两个相同，但速度要快得多，如何？

---

我在我的系统上尝试了 Z boson 的程序以及我的 2 个实现，得到了以下结果：

> g++ -O2 -fopenmp main.cpp
> ./a.out
mul         time 1.33, 18.08 GB/s
mul_SSE     time 1.13, 21.24 GB/s
mul_SSE_NT  time 1.51, 15.88 GB/s
mul_SSE_OMP time 0.79, 30.28 GB/s
mul_SSE_v2  time 1.12, 21.49 GB/s
mul_v2      time 1.26, 18.99 GB/s
mul_asm     time 1.12, 21.50 GB/s
mul_asm2    time 1.09, 22.08 GB/s

Answer 1

以前的基准测试有 a major bug in the timing function I used。这严重低估了没有矢量化和其他测量的带宽。此外，还有另一个问题是高估了阵列上已读取但未写入的带宽due to COW。最后，我使用的最大带宽不正确。我已通过更正更新了我的答案，并在此答案的末尾留下了旧答案。

---

您的操作受内存带宽限制。这意味着 CPU 大部分时间都在等待缓慢的内存读取和写入。可以在这里找到一个很好的解释：Why vectorizing the loop does not have performance improvement。

但是，我不得不稍微不同意那个答案中的一个陈述。

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

事实上，即使在内存带宽受限的操作中，矢量化、展开和多线程也可以显着增加带宽。原因是难以获得最大的内存带宽。可以在这里找到一个很好的解释：https://stackoverflow.com/a/25187492/2542702。

我的其余答案将展示矢量化和多线程如何接近最大内存带宽。

我的测试系统：Ubuntu 16.10，Skylake (i7-6700HQ@2.60GHz)，32GB RAM，双通道 DDR4@2400 GHz。我系统的最大带宽为 38.4 GB/s。

根据下面的代码，我生成了下表。我使用 OMP_NUM_THREADS 例如设置线程数export OMP_NUM_THREADS=4。效率为bandwidth/max_bandwidth。

-O2 -march=native -fopenmp
Threads Efficiency
1       59.2%
2       76.6%
4       74.3%
8       70.7%

-O2 -march=native -fopenmp -funroll-loops
1       55.8%
2       76.5%
4       72.1%
8       72.2%

-O3 -march=native -fopenmp
1       63.9%
2       74.6%
4       63.9%
8       63.2%

-O3 -march=native -fopenmp -mprefer-avx128
1       67.8%
2       76.0%
4       63.9%
8       63.2%

-O3 -march=native -fopenmp -mprefer-avx128 -funroll-loops
1       68.8%
2       73.9%
4       69.0%
8       66.8%

由于测量中存在不确定性，经过多次迭代运行，我得出以下结论：

• 单线程标量操作获得超过 50% 的带宽。
• 两个线程标量操作获得最高带宽。
• 单线程向量操作比单线程标量操作更快。
• 单线程 SSE 操作比单线程 AVX 操作更快。
• 展开没有帮助。
• 展开单线程操作比不展开要慢。
• 线程数多于内核数（超线程）会降低带宽。

提供最佳带宽的解决方案是具有两个线程的标量操作。

我用来基准测试的代码：

#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <omp.h>

#define N 10000000
#define R 100

void mul(double *a, double *b) {
  #pragma omp parallel for
  for (int i = 0; i<N; i++) a[i] *= b[i];
}

int main() {
  double maxbw = 2.4*2*8; // 2.4GHz * 2-channels * 64-bits * 1-byte/8-bits 
  double mem = 3*sizeof(double)*N*R*1E-9; // GB

  double *a = (double*)malloc(sizeof *a * N);
  double *b = (double*)malloc(sizeof *b * N);

  //due to copy-on-write b must be initialized to get the correct bandwidth
  //also, GCC will convert malloc + memset(0) to calloc so use memset(1)
  memset(b, 1, sizeof *b * N);

  double dtime = -omp_get_wtime();
  for(int i=0; i<R; i++) mul(a,b);
  dtime += omp_get_wtime();
  printf("%.2f s, %.1f GB/s, %.1f%%\n", dtime, mem/dtime, 100*mem/dtime/maxbw);

  free(a), free(b);
}

---

有计时错误的旧解决方案

内联汇编的现代解决方案是使用内部函数。仍然存在需要内联汇编的情况，但这不是其中之一。

内联汇编方法的一个内在解决方案很简单：

void mul_SSE(double*  a, double*  b) {
  for (int i = 0; i<N/2; i++) 
      _mm_store_pd(&a[2*i], _mm_mul_pd(_mm_load_pd(&a[2*i]),_mm_load_pd(&b[2*i])));
}

让我定义一些测试代码

#include <x86intrin.h>
#include <string.h>
#include <stdio.h>
#include <x86intrin.h>
#include <omp.h>

#define N 1000000
#define R 1000

typedef __attribute__(( aligned(32)))  double aligned_double;
void  (*fp)(aligned_double *a, aligned_double *b);

void mul(aligned_double* __restrict a, aligned_double* __restrict b) {
  for (int i = 0; i<N; i++) a[i] *= b[i];
}

void mul_SSE(double*  a, double*  b) {
  for (int i = 0; i<N/2; i++) _mm_store_pd(&a[2*i], _mm_mul_pd(_mm_load_pd(&a[2*i]),_mm_load_pd(&b[2*i])));
}

void mul_SSE_NT(double*  a, double*  b) {
  for (int i = 0; i<N/2; i++) _mm_stream_pd(&a[2*i], _mm_mul_pd(_mm_load_pd(&a[2*i]),_mm_load_pd(&b[2*i])));
}

void mul_SSE_OMP(double*  a, double*  b) {
  #pragma omp parallel for
  for (int i = 0; i<N; i++) a[i] *= b[i];
}

void test(aligned_double *a, aligned_double *b, const char *name) {
  double dtime;
  const double mem = 3*sizeof(double)*N*R/1024/1024/1024;
  const double maxbw = 34.1;
  dtime = -omp_get_wtime();
  for(int i=0; i<R; i++) fp(a,b);
  dtime += omp_get_wtime();
  printf("%s \t time %.2f s, %.1f GB/s, efficency %.1f%%\n", name, dtime, mem/dtime, 100*mem/dtime/maxbw);
}

int main() {
  double *a = (double*)_mm_malloc(sizeof *a * N, 32);
  double *b = (double*)_mm_malloc(sizeof *b * N, 32);

  //b must be initialized to get the correct bandwidth!!!
  memset(a, 1, sizeof *a * N);
  memset(b, 1, sizeof *a * N);

  fp = mul,         test(a,b, "mul        ");
  fp = mul_SSE,     test(a,b, "mul_SSE    ");
  fp = mul_SSE_NT,  test(a,b, "mul_SSE_NT ");
  fp = mul_SSE_OMP, test(a,b, "mul_SSE_OMP");

  _mm_free(a), _mm_free(b);
}

现在是第一个测试

g++ -O2 -fopenmp test.cpp
./a.out
mul              time 1.67 s, 13.1 GB/s, efficiency 38.5%
mul_SSE          time 1.00 s, 21.9 GB/s, efficiency 64.3%
mul_SSE_NT       time 1.05 s, 20.9 GB/s, efficiency 61.4%
mul_SSE_OMP      time 0.74 s, 29.7 GB/s, efficiency 87.0%

所以-O2 不向量化循环，我们看到内在 SSE 版本比普通 C 解决方案 mul 快得多。 efficiency = bandwith_measured/max_bandwidth 我的系统的最大值为 34.1 GB/s。

第二次测试

g++ -O3 -fopenmp test.cpp
./a.out
mul              time 1.05 s, 20.9 GB/s, efficiency 61.2%
mul_SSE          time 0.99 s, 22.3 GB/s, efficiency 65.3%
mul_SSE_NT       time 1.01 s, 21.7 GB/s, efficiency 63.7%
mul_SSE_OMP      time 0.68 s, 32.5 GB/s, efficiency 95.2%

-O3 对循环进行矢量化处理，而内部函数基本上没有任何优势。

第三次测试

g++ -O3 -fopenmp -funroll-loops test.cpp
./a.out
mul              time 0.85 s, 25.9 GB/s, efficency 76.1%
mul_SSE          time 0.84 s, 26.2 GB/s, efficency 76.7%
mul_SSE_NT       time 1.06 s, 20.8 GB/s, efficency 61.0%
mul_SSE_OMP      time 0.76 s, 29.0 GB/s, efficency 85.0%

使用-funroll-loops GCC 将循环展开八次，除了非临时存储解决方案和 OpenMP 解决方案没有真正优势外，我们看到了显着的改进。

在展开循环之前，mul 和 -O3 的程序集是

    xor     eax, eax
.L2:
    movupd  xmm0, XMMWORD PTR [rsi+rax]
    mulpd   xmm0, XMMWORD PTR [rdi+rax]
    movaps  XMMWORD PTR [rdi+rax], xmm0
    add     rax, 16
    cmp     rax, 8000000
    jne     .L2
    rep ret

对于-O3 -funroll-loops，mul 的程序集是：

   xor     eax, eax
.L2:
    movupd  xmm0, XMMWORD PTR [rsi+rax]
    movupd  xmm1, XMMWORD PTR [rsi+16+rax]
    mulpd   xmm0, XMMWORD PTR [rdi+rax]
    movupd  xmm2, XMMWORD PTR [rsi+32+rax]
    mulpd   xmm1, XMMWORD PTR [rdi+16+rax]
    movupd  xmm3, XMMWORD PTR [rsi+48+rax]
    mulpd   xmm2, XMMWORD PTR [rdi+32+rax]
    movupd  xmm4, XMMWORD PTR [rsi+64+rax]
    mulpd   xmm3, XMMWORD PTR [rdi+48+rax]
    movupd  xmm5, XMMWORD PTR [rsi+80+rax]
    mulpd   xmm4, XMMWORD PTR [rdi+64+rax]
    movupd  xmm6, XMMWORD PTR [rsi+96+rax]
    mulpd   xmm5, XMMWORD PTR [rdi+80+rax]
    movupd  xmm7, XMMWORD PTR [rsi+112+rax]
    mulpd   xmm6, XMMWORD PTR [rdi+96+rax]
    movaps  XMMWORD PTR [rdi+rax], xmm0
    mulpd   xmm7, XMMWORD PTR [rdi+112+rax]
    movaps  XMMWORD PTR [rdi+16+rax], xmm1
    movaps  XMMWORD PTR [rdi+32+rax], xmm2
    movaps  XMMWORD PTR [rdi+48+rax], xmm3
    movaps  XMMWORD PTR [rdi+64+rax], xmm4
    movaps  XMMWORD PTR [rdi+80+rax], xmm5
    movaps  XMMWORD PTR [rdi+96+rax], xmm6
    movaps  XMMWORD PTR [rdi+112+rax], xmm7
    sub     rax, -128
    cmp     rax, 8000000
    jne     .L2
    rep ret

第四次测试

g++ -O3 -fopenmp -mavx test.cpp
./a.out
mul              time 0.87 s, 25.3 GB/s, efficiency 74.3%
mul_SSE          time 0.88 s, 24.9 GB/s, efficiency 73.0%
mul_SSE_NT       time 1.07 s, 20.6 GB/s, efficiency 60.5%
mul_SSE_OMP      time 0.76 s, 29.0 GB/s, efficiency 85.2%

现在非内在函数是最快的（不包括 OpenMP 版本）。

因此在这种情况下没有理由使用内在函数或内联汇编，因为我们可以通过适当的编译器选项（例如-O3、-funroll-loops、-mavx）获得最佳性能。

测试系统：Ubuntu 16.10，Skylake (i7-6700HQ@2.60GHz)，32GB RAM。最大内存带宽（34.1 GB/s）https://ark.intel.com/products/88967/Intel-Core-i7-6700HQ-Processor-6M-Cache-up-to-3_50-GHz

---

这是另一个值得考虑的解决方案。 The cmp instruction is not necessary 如果我们从 -N 数到零并以N+i 访问数组。 GCC 早就应该解决这个问题了。它消除了一条指令（尽管由于宏操作融合， cmp 和 jmp 通常算作一个微操作）。

void mul_SSE_v2(double*  a, double*  b) {
  for (ptrdiff_t i = -N; i<0; i+=2)
    _mm_store_pd(&a[N + i], _mm_mul_pd(_mm_load_pd(&a[N + i]),_mm_load_pd(&b[N + i])));

与-O3 组装

mul_SSE_v2(double*, double*):
    mov     rax, -1000000
.L9:
    movapd  xmm0, XMMWORD PTR [rdi+8000000+rax*8]
    mulpd   xmm0, XMMWORD PTR [rsi+8000000+rax*8]
    movaps  XMMWORD PTR [rdi+8000000+rax*8], xmm0
    add     rax, 2
    jne     .L9
    rep ret
}

这种优化只可能有助于数组适合，例如L1 缓存，即不从主存读取。

---

我终于找到了一种方法来获得不生成 cmp 指令的普通 C 解决方案。

void mul_v2(aligned_double* __restrict a, aligned_double* __restrict b) {
  for (int i = -N; i<0; i++) a[i] *= b[i];
}

然后从像 mul_v2(&a[N],&b[N]) 这样的单独对象文件中调用函数，所以这可能是最好的解决方案。但是，如果您从与在 GCC 中定义的对象文件（翻译单元）相同的对象文件（翻译单元）中调用该函数，则会再次生成 cmp 指令。

还有，

void mul_v3(aligned_double* __restrict a, aligned_double* __restrict b) {
  for (int i = -N; i<0; i++) a[N+i] *= b[N+i];
}

仍会生成cmp 指令并生成与mul 函数相同的程序集。

---

函数mul_SSE_NT 很傻。它使用非临时存储，仅在仅写入内存时才有用，但由于函数读取和写入同一地址，非临时存储不仅无用，而且结果较差。

---

此答案的先前版本获得了错误的带宽。原因是数组没有初始化。

Answer 2

我想为这个问题添加另一个观点。 如果没有内存限制，SIMD 指令会大大提高性能。但是在当前示例中，内存加载和存储操作太多，CPU 计算太少。因此 CPU 可以及时处理传入的数据，而无需使用 SIMD。 如果您使用其他类型的数据（例如 32 位浮点数）或更复杂的算法，内存吞吐量不会限制 CPU 性能，使用 SIMD 会带来更多优势。

Answer 3

你的 asm 代码真的没问题。 不是你衡量它的方式。 正如我在 cmets 中指出的那样，您应该：

a) 使用更多的迭代次数 - 100 万次对于现代 CPU 来说不算什么

b) 使用 HPT 进行测量

c) 使用 RDTSC 或 RDTSCP 计算实际 CPU 时钟数

另外你为什么害怕-O3 opt？不要忘记为您的平台构建代码，因此请使用 -march=native。如果您的 CPU 支持 AVX 或 AVX2 编译器将借此机会生成更好的代码。

接下来 - 如果你知道你的代码，给编译器一些关于别名和对齐的提示。

这是我的 mul_c 版本 - 是的，它是 GCC 特定的，但你表明你使用了 GCC

void mul_c(double* restrict a, double* restrict b)
{
   a = __builtin_assume_aligned (a, 16);
   b = __builtin_assume_aligned (b, 16);

    for (int i = 0; i != 1000000; ++i)
    {
        a[i] = a[i] * b[i];
    }
}

它将产生：

mul_c(double*, double*):
        xor     eax, eax
.L2:
        movapd  xmm0, XMMWORD PTR [rdi+rax]
        mulpd   xmm0, XMMWORD PTR [rsi+rax]
        movaps  XMMWORD PTR [rdi+rax], xmm0
        add     rax, 16
        cmp     rax, 8000000
        jne     .L2
        rep ret

如果你有 AVX2 并确保数据是 32 字节对齐的，它将变成

mul_c(double*, double*):
        xor     eax, eax
.L2:
        vmovapd ymm0, YMMWORD PTR [rdi+rax]
        vmulpd  ymm0, ymm0, YMMWORD PTR [rsi+rax]
        vmovapd YMMWORD PTR [rdi+rax], ymm0
        add     rax, 32
        cmp     rax, 8000000
        jne     .L2
        vzeroupper
        ret

因此，如果编译器可以为您完成，则无需手工制作 asm ;)