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

Question

使用 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的单核上执行指令。我的酷睿 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

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

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

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

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

但有时，即使是最好的编译器也无法给出令人满意的结果(始终检查您的程序集的性能关键部分)。

他们需要帮助，有时您需要自己编写汇编。

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

首先，您可以实现经典的优化技巧:

1. 常量和别名

通过 __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彼此。

1. 对齐信息:

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

#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(加载未对齐)。

1. 展开内部循环(如果可以的话):

如果您知道您的问题大小是 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
...

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

1. 内部函数

如果你想确保你有正确的汇编，那么你必须去内在函数/汇编。

一些简单的东西，比如:

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

1. 非临时存储:作为最后的优化，您可以使用 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));
    }
}

这段代码在我的至强 E5-1620 v3 上达到了 34 GB/s，DDR4 内存 @ 2133 MHz，而开始时的简单解决方案是 33 GB/s。

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

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