What's the fastest way to multi-threaded SIMD operation?

If the problem is big, you should be multithreaded.

You can choose openmp or pthread, they will give you similar levels of performance (maybe slightly better with pthread, but this will not be portable and more difficult to maintain).

Your code will be bandwidth limited, not computed.

To achieve maximum throughput, you need to interleave independent memory operations using multithreading.

very simple solution like

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

      

        
        
        
      

    

will probably give you acceptable performance on all systems with each compiler.

In fact, letting the compiler optimize will likely give you good results and will most likely help you write readable code.

But sometimes even the best compiler does not give satisfactory results (always check your build on performance-critical sections).

They need help and sometimes you need to write the assembly yourself.

Here is the path I will follow to optimize this loop until I get the results I want.

First, there are the classic optimization tricks you can implement:

• constant and alias

Provide a constant and prevent aliasing with the __restrict keyword:

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

      

        
        
        
      

    

This will help the compiler since it will know that a, b and c cannot alias each other.

2. alignment information :

Tell the compiler that your pointers are correctly aligned

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

      

        
        
        
      

    

This will also help the compiler generate vload instead of vloadu (load unaligned).

3. Expand inner loops (if possible):

If you know your problem size, if multiple of 256 bits, you can even expand the inner loop:

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

      

        
        
        
      

    

with this code, clang 4.0 gives a pretty neat build:

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

      

        
        
        
      

    

For some reason, you need to tweak your attributes and pragmas to have the same result with other compilers.

4. Intrinsics

If you want to make sure you have the correct assembly, you should go to intrinsics / assembly.

Something simple:

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

      

        
        
        
      

    

5. Non-temporary storage : As a final optimization, you can use the non-temporary hint in the storage statement, since another iteration of the loop will not read the value you just wrote:

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

      

        
        
        
      

    

which gives you this assembly:

.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

      

        
        
        
      

    

If you are worried about the cmp instruction every step of the way, you can unwrap more steps in your loop, but branch prediction does a pretty good job on modern CPUs

[EDIT: add pthread] As stated above, pthread is a bit painful to manage ... Here is a complete functional example with 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));
    }
}

      

        
        
        
      

    

This code goes up to 34GB / s on my Xeon E5-1620 v3 with DDR4 at 2133MHz and the simple solution at the start is 33GB / s.

All these efforts to save 3% :). But sometimes this 3% can be critical.

Note that memory initialization must be done by the same kernel that will perform the computation (especially for NUMA systems) to avoid page migrations.