显式多线程 SIMD 操作的最快方法是什么？答案

【问题标题】：What is the fastest way for a multithread SIMD operation explicitly?显式多线程 SIMD 操作的最快方法是什么？
【发布时间】：2017-05-25 21:50:10
【问题描述】：

使用intrinsics 是SIMDizing 的常用方法。例如，我可以通过_mm256_add_epi32 对八个整数执行一条加法指令。添加后需要两个_mm256_load_si256和一个_mm256_store_si256，如下：

__m256i vec1 = _mm256_load_si256((__m256i *)&A[0]); // almost 5 cycles
__m256i vec2 = _mm256_load_si256((__m256i *)&B[0]); // almost 5 cycles
__m256i vec3 = _mm256_add_epi32( vec1 , vec2); // almost 1 cycle
_mm256_store_si256((__m256i *)&C[0], vec3); // almost 5

它在 CPU 的单核上执行指令。我的 Core i7 有 8 个核心（4 个实心）；我想像这样将操作发送到所有核心：

int i_0, i_1, i_2, i_3, i_4, i_5, i_6, i_7 ; // These specify the values in memory
//core 0
__m256i vec1_0 = _mm256_load_si256((__m256i *)&A[i_0]);  
__m256i vec2_0 = _mm256_load_si256((__m256i *)&B[i_0]); 
__m256i vec3_0 = _mm256_add_epi32( vec1 , vec2); 
_mm256_store_si256((__m256i *)&C[i_0], vec3_0);

//core 1
__m256i vec1_1 = _mm256_load_si256((__m256i *)&A[i_1]);
__m256i vec2_1 = _mm256_load_si256((__m256i *)&B[i_1]);
__m256i vec3_1 = _mm256_add_epi32( vec1 , vec2);
_mm256_store_si256((__m256i *)&C[i_1], vec3_1);

//core 2
__m256i vec1_2 = _mm256_load_si256((__m256i *)&A[i_2]);
__m256i vec2_2 = _mm256_load_si256((__m256i *)&B[i_2]);
__m256i vec3_2 = _mm256_add_epi32( vec1 , vec2);
_mm256_store_si256((__m256i *)&C[i_2], vec3_2);

//core 3
__m256i vec1_3 = _mm256_load_si256((__m256i *)&A[i_3]);
__m256i vec2_3 = _mm256_load_si256((__m256i *)&B[i_3]);
__m256i vec3_3 = _mm256_add_epi32( vec1 , vec2);
_mm256_store_si256((__m256i *)&C[i_3], vec3_3);

//core 4
__m256i vec1_4 = _mm256_load_si256((__m256i *)&A[i_4]);
__m256i vec2_4 = _mm256_load_si256((__m256i *)&B[i_4]);
__m256i vec3_4 = _mm256_add_epi32( vec1 , vec2);
_mm256_store_si256((__m256i *)&C[i_4], vec3_4);

//core 5
__m256i vec1_5 = _mm256_load_si256((__m256i *)&A[i_5]);
__m256i vec2_5 = _mm256_load_si256((__m256i *)&B[i_5]);
__m256i vec3_5 = _mm256_add_epi32( vec1 , vec2);
_mm256_store_si256((__m256i *)&C[i_5, vec3_5);

//core 6
__m256i vec1_6 = _mm256_load_si256((__m256i *)&A[i_6]);
__m256i vec2_6 = _mm256_load_si256((__m256i *)&B[i_6]);
__m256i vec3_6 = _mm256_add_epi32( vec1 , vec2);
_mm256_store_si256((__m256i *)&C[i_6], vec3_6);

//core 7
__m256i vec1_7 = _mm256_load_si256((__m256i *)&A[i_7]);
__m256i vec2_7 = _mm256_load_si256((__m256i *)&B[i_7]);
__m256i vec3_7 = _mm256_add_epi32( vec1 , vec2);
_mm256_store_si256((__m256i *)&C[i_7], vec3_7);

POSIX Thread 可用，openMP 在这种情况下也很有用。但是，与此操作的几乎5+5+1 cyles 相比，创建和维护线程花费了太多时间。因为，所有数据都是依赖的，所以我不需要查看共享内存。实现此操作的最快显式方法是什么？

因此，我从事 GPP 工作，GPU 可能不是答案。我还想实现一个库，因此编译器基础解决方案可能是一个挑战者。这个问题对于多线程来说已经足够大了。这是为了我的研究，因此我可以改变问题以适应这个概念。我想实现一个库并将其与 OpenMP 等其他解决方案进行比较，希望我的库将比其他当前解决方案更快。 GCC 6.3/clang 3.8, Linux Mint, Skylake

提前致谢。

【问题讨论】：

除非这是在一个深度循环的底部，在另一个循环被多次调用，并且你有一个显示它是瓶颈的配置文件，那么你可以选择的任何方法都将是最快的。在不知道是否确实存在问题的情况下编写最快的代码时请小心。如果这是需要更快的代码，那么继续进行实验，看看分析器说了什么。我什至不知道上面的代码会按原样使用多个 CPU - 相反，我相信它会将它们全部排在同一个核心上。
@MichaelDorgan，谢谢，我把问题改得更笼统了。它不是在一个循环中，但是，它可以是。我正在为我的应用程序实现一个多线程 SIMD 库，它是我的问题的简化版本。
一个问题是所有内核都在竞争相同的内存和 L3 和/或 L4 缓存。如果进程的内存带宽受限，只有 1 或 2 个内核，那么使用额外的内核将无济于事。
把你的问题分解成块。如果在将其拆分为多个内核时，您的内存访问模式是分散的，那么您从多线程中获得的好处很少。相反，如果您在内存中进行了严格的操作，那么将其划分到多个内核上并在较小的块上操作可能会有所帮助。您可能需要添加一些预取指令来帮助提前准备加载。
只有当问题足够大时，在多核上拆分计算才有意义。在你的情况下，它绝对没有。所以你需要使用操作系统提供的线程功能。如果您想减少线程创建开销，请考虑使用线程池。

标签： c multithreading x86 simd intrinsics

【解决方案1】：

如果你的问题很大，你必须多线程。

您可以选择 openmp 或 pthread，它们会为您提供相似的性能水平（使用 pthread 可能会更好一些，但这将是不可移植的并且维护起来更复杂）。

您的代码将受带宽限制，绝对不受计算限制。

为了达到最大吞吐量，需要通过多线程来交错独立的内存操作。

一个非常简单的解决方案，比如

extern "C" void add(int* a, int* b, int* c, int N) {
    #pragma omp parallel for
    for(int i = 0; i < N; ++i) {
        a[i] = b[i] + c[i];
    }
}

可能会在所有系统上使用每个编译器为您提供可接受的性能。

事实上，让编译器优化可能会给你带来很好的性能，并且肯定会帮助你编写可读的代码。

但有时，即使是最好的编译器也不能给出令人满意的结果（总是检查你的程序集的性能关键部分）。

他们需要帮助，有时您需要自己编写程序集。

这是我将遵循的优化此循环的路径，直到获得我想要的结果。

首先，您可以实施一些经典的优化技巧：

常量和别名

通过 __restrict 关键字提供常量并防止别名：

extern "C" void add(int* __restrict a, const int* __restrict b, const int* __restrict c, int N) {
    #pragma omp parallel for
    for(int i = 0; i < N; ++i) {
        a[i] = b[i] + c[i];
    }
}

这将有助于编译器，因为它会知道 a、b 和 c 不能alias 彼此。

对齐信息：

告诉编译器你的指针正确对齐

#define RESTRICT __restrict

    typedef __attribute__((aligned(32))) int* intptr;

    extern "C" void add(intptr RESTRICT a, const intptr RESTRICT b, const intptr RESTRICT c, int N) {
        #pragma omp parallel for
        for(int i = 0; i < N; ++i) {
            a[i] = b[i] + c[i];
        }
    }

这也将有助于编译器生成 vload 指令而不是 vloadu（未对齐加载）。

展开内部循环（如果可以的话）：

如果你知道你的问题大小是 256 位的倍数，你甚至可以展开一个内部循环：

#define RESTRICT __restrict

typedef __attribute__((aligned(32))) int* intptr;

extern "C" void add(intptr RESTRICT a, const intptr RESTRICT b, const intptr RESTRICT c, int N) {
    #pragma omp parallel for
    for(int i = 0; i < N; i += 8) {
        #pragma unroll
        for(int k = 0; k < 8; ++k)
        a[i+k] = b[i+k] + c[i+k];
    }
}

使用该代码，clang 4.0 提供了相当简洁的组装：

...
 vmovdqu ymm0, ymmword ptr [rdx + 4*rcx]
 vpaddd  ymm0, ymm0, ymmword ptr [rsi + 4*rcx]
 vmovdqu ymmword ptr [rdi + 4*rcx], ymm0
...

由于某些原因，您需要调整您的属性和编译指示以与其他编译器获得相同的结果。

内在函数

如果你想确定你有正确的程序集，那么你必须去内在/程序集。

一些简单的东西：

#define RESTRICT __restrict

typedef __attribute__((aligned(32))) int* intptr;

extern "C" void add(intptr RESTRICT a, const intptr RESTRICT b, const intptr RESTRICT c, int N) {
    #pragma omp parallel for
    for(int i = 0; i < N; i += 8) {
        __m256i vb = _mm256_load_si256((__m256i*) (b + i));
        __m256i vc = _mm256_load_si256((__m256i*) (c + i));
        _mm256_store_si256((__m256i*) (a + i), _mm256_add_epi32(vb, vc));
    }
}

非临时存储：作为最后的优化，您可以在 store 指令上使用 non-temporal hint，因为循环的其他迭代不会读取您刚刚写入的值：

typedef __attribute__((aligned(32))) int* intptr;
extern "C" void add(intptr RESTRICT a, const intptr RESTRICT b, const intptr RESTRICT c, int N) {
    #pragma omp parallel for
    for(int i = 0; i < N; i += 8) {
        __m256i vb = _mm256_load_si256((__m256i*) (b + i));
        __m256i vc = _mm256_load_si256((__m256i*) (c + i));
        _mm256_stream_si256((__m256i*) (a + i), _mm256_add_epi32(vb, vc));
    }
}

它为您提供了该程序集：

.L3:
        vmovdqa ymm0, YMMWORD PTR [rdx+rax]
        vpaddd  ymm0, ymm0, YMMWORD PTR [rsi+rax]
        vmovntdq        YMMWORD PTR [rdi+rax], ymm0
        add     rax, 32
        cmp     rcx, rax
    jne     .L3
    vzeroupper

如果您对每一步的 cmp 指令感到担心，您可以在循环中展开更多步骤，但 branch prediction 在现代处理器上做得相当不错

[编辑：添加 pthread] 如上所述，pthread 管理起来有点痛苦...... 这是一个完整的 pthread 示例：

#include <pthread.h>
#include <cstdlib>
#include <cstdio>
#include <immintrin.h>

typedef struct AddStruct {
    int *a, *b, *c;
    int N;
} AddStruct_t;

void* add(void* s);

int main() {
    const int N = 1024*1024*32; // out of cache
    int *a, *b, *c;
    int err;
    err = posix_memalign((void**) &a, 32, N*sizeof(int));
    err = posix_memalign((void**) &b, 32, N*sizeof(int));
    err = posix_memalign((void**) &c, 32, N*sizeof(int));
    for(int i = 0; i < N; ++i) {
        a[i] = 0;
        b[i] = 1;
        c[i] = i;
    }
int slice = N / 8;
pthread_t threads[8];
AddStruct_t arguments[8];
for(int i = 0; i < 8; ++i) {
    arguments[i].a = a + slice * i;
    arguments[i].b = b + slice * i;
    arguments[i].c = c + slice * i;
    arguments[i].N = slice;
}

for(int i = 0; i < 8; ++i) {
    if(pthread_create(&threads[i], NULL, add, &arguments[i])) {
        fprintf(stderr, "ERROR CREATING THREAD %d\n", i);
        abort();
    }
   }

for(int i = 0; i < 8; ++i) {
    pthread_join(threads[i], NULL);
}

for(int i = 0; i < N; ++i) {
    if(a[i] != i + 1) {
        fprintf(stderr, "ERROR AT %d: expected %d, actual %d\n", i, i+1, a[i]);
        abort();
    }
}

fprintf(stdout, "OK\n");
}

void* add(void* v) {
    AddStruct_t* s = (AddStruct_t*) v;
    for(int i = 0; i < s->N; i += 8) {
        __m256i vb = _mm256_load_si256((__m256i*) (s->b + i));
        __m256i vc = _mm256_load_si256((__m256i*) (s->c + i));
        _mm256_stream_si256((__m256i*) (s->a + i), _mm256_add_epi32(vb, vc));
    }
}

这段代码在我的 Xeon E5-1620 v3 上实现了 34 GB/s，DDR4 内存 @ 2133 MHz，而一开始的简单解决方案是 33 GB/S。

所有这些努力都是为了节省 3% :)。但有时这 3% 可能很关键。

请注意，内存初始化应由执行计算的同一内核执行（尤其是 NUMA 系统）以避免页面迁移。

【讨论】：

谢谢，一个很好的答案。你也可以添加一些显式的多线程吗？可移植性不是我关心的问题，因为我在 Skylake 上执行代码。我使用了pthreads，但我的解决方案并不好。事实上，我可能会更改编译器并希望库独立工作。
我添加了一个带有 pthread 的示例，但是，为了代码的可读性，我真的会考虑 openmp。不管怎样，做你的性能测量，然后选择你喜欢的多线程库
您在答案开头的简单解决方案就是您所需要的。对于小 N，您不想使用多线程，而对于大 N，它的内存带宽受限，因此展开、对齐、限制实际上是无用的。
顺便说一句，您可以使用for simd aligned(a,b,c:32)（我的语法可能不正确）进行对齐。
GCC 的 -funroll-loops 在使用 OpenMP 的情况下比没有使用 OpenMP 时优化得更差（我通过将函数移动到单独的目标文件中发现了这一点，并且它变得更快）。展开是您仍然必须使用 GCC 手动完成的操作才能获得良好的结果。同样，尽管这些优化仅在操作不受内存带宽限制的情况下才有用。 OP的问题太宽泛了。

【解决方案2】：

实现最快的方法是：

void add_ints(int *vec1, int *vec2, int *vec3 int n){
 int i; 
#pragma simd
for (i=0; i<n; i++){
  vec3[i] = vec1[i] + vec2[i] ;
}

“自己动手”是否更快值得调查。但是“自己动手”可能更容易出错......这使得实施速度变慢。

对于这些简单的问题，人们会期望编译器编写者足够老练，能够理解简单问题的最快解决方案，而且通常他们甚至可以很好地找到复杂问题的最快解决方案......而且#pragma 的使用有助于他们。

其次；与单核上的直接“SIMD”相比，我很少发现“SIMD 并行”在 IO 驱动的问题（例如 ^this^）上工作得更快。
我通常会实现略低于 1600 MB/秒的吞吐量，这在 1600 内存上似乎相当不错。
除非 GPU 的 IO 带宽高于 1600 MB/秒，否则在单个主机内核上可能会更好，并在需要更多数学/IO 时使用 GPU。

但是，您可以而且应该亲自尝试一下。（是的...以下示例不在 icc 网站上）

#pragma omp parallel for simd schedule(static,10) {
  for (i=0; i<N; i++) { vec3[i] = vec1[i] + vec2[i]; }
}

在您掌握了简单的方法之后，您可以测量“您拥有的滚动”在使用单核和多核的 -O3 编译器上的性能有多好。

向量的另一个选择是 CILK+ 。当具有 MATLAB 或 Fortran 背景的人尤其如此，因为向量和矩阵/数组结构非常相似。

基本上，SIMD 内在函数在早期就“流行”了，一旦编译器和 OpenMP 将它们引入到它们的内部，那么仅在编译器无法提供向量机的情况下使用内在函数似乎会更好- 为您编码。

【讨论】：

我所知道的唯一支持#pragma sims 的编译器是ICC，这不是问这个问题的人正在使用或希望针对的目标。它也应该是完全没有必要的。如果您使用适当的编译器开关，优化器将做出正确的决定，而不需要分散在整个代码库中的不可移植的#pragmas。
#pragma omp simd 从 ~4.9-5.1 版本开始受到 GCC 的支持。所以它被 icc、ifort 和 gcc 以及一些基于 llvm 的支持。我知道的唯一一个完全不支持的 x86 编译器是 MS 编译器。
#pragmas 很容易被编译器忽略，可能有 Gcc 版本的开关？是否甚至需要编译指示都需要查看矢量报告。 OP 以“最容易实现”结束了这个问题，即使我能做到，编译指示也很容易。但是如果编译器优化了代码，那么就不需要编译指示。执行汇编程序类型加载似乎更难实现。 OP 有一个寻找问题的解决方案，并且应该根据可靠的事实将自定义库放在一起。因此尝试多种方法是合理的。
谢谢，您的解决方案很棒，但我想明确地实现它。
这是有道理的@FackedDeveloper 理想情况下，您可以回复它的运行速度。