为什么这种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"
    );
}

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

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

mul_c: 1 ms
mul_asm: 1 ms

[a lot of doubles...]

我预计SIMD实现的速度至少是乘法/内存指令数量的一半(0 ms)的两倍。

因此，我的问题是：为什么不能比普通的C/C++实现更快，而SIMD实现只执行乘法/内存指令量的一半？

以下是完整的程序：

#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"
    );
}

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

更新

正如答案/评论所建议的那样，我实现了另一种测量执行时间的方法：

#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玻色子的程序以及我的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

您的asm代码真的很好。不是你测量它的方式。正如我在评论中指出的那样，你应该：

( a)使用更多的迭代方式--100万对于现代CPU来说是没有意义的。

( b)使用HPT进行测量

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

另外，你为什么害怕-O3选项？不要忘记为您的平台构建代码，所以使用-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 ;)

Answer 2

我还想对这个问题补充另一点看法。如果没有内存绑定限制，SIMD指令会大大提高性能。但是在当前的例子中，内存加载和存储操作太多，CPU计算太少。因此，CPU可以在不使用SIMD的情况下及时处理传入的数据。如果您使用另一种类型的数据(例如32位浮点数)或更复杂的算法，内存吞吐量不会限制CPU性能，使用SIMD将提供更多的优势。