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

Question

使用intrinsics是SIMD化的常用方法。例如，我可以通过 _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的单核上执行指令。我的酷睿 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 线程可用，openMP 在这种情况下也很有用。但是，与此操作的几乎 5+5+1 个周期相比，创建和维护线程需要花费太多时间。因为，所有数据都是相关的，所以我不需要看共享内存。实现此操作最快的显式方法是什么？

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

提前致谢。

Answer 1

实施最快的方法是：

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] ;
}

"roll your own" 是否更快值得研究。但是 "rolling your own" 更容易出错......这使得实施起来更慢。

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

其次；我很少发现 'SIMD parallel' 在 IO 驱动问题上工作得更快的情况，例如 ^this^，与单核上的直接 'SIMD' 相比。
我经常达到略低于 1600 MB/second 的吞吐量，这在 1600 内存上似乎相当不错。
除非 GPU 的 IO 带宽高于 1600 MB/sec，否则您在单个主机内核上可能会更好，并在需要更多 math/IO 时使用 GPU。

但是您可以而且应该亲自尝试一下。（是的...下面的例子来自icc网站）

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

掌握简单方法后，您可以测量 "roll you own" 比使用单核和多核的 -O3 编译器的性能好多少。

向量的另一个选择是 CILK+。当一个人来自 MATLAB 或 Fortran 背景时尤其如此，因为向量和 matrix/array 构造非常相似。

基本上 SIMD 内在函数是 'in vogue' 早期的，一旦编译器和 OpenMP 将它们纳入其内部，那么仅在编译器无法使用的情况下保留使用内在函数似乎更好为您提供矢量机器码。

Answer 2

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

您可以选择 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));
    }
}

非临时存储：作为最后的优化，您可以在存储指令上使用 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 和 DDR4 内存 @ 2133 MHz 上实现了 34 GB/s，而开始时的简单解决方案是 33 GB/S。

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

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

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

What is the fastest way for a multithread SIMD operation explicitly?

c

x86

multithreading

simd

intrinsics