Правильный способ суммировать два массива с SSE2 SIMD в C++

Question

Давайте начнем в том числе следующие:

#include <vector> 
#include <random> 
using namespace std; 

Теперь предположим, что один имеет следующие три std:vector<float>:

N = 1048576; 
vector<float> a(N); 
vector<float> b(N); 
vector<float> c(N); 

default_random_engine randomGenerator(time(0)); 
uniform_real_distribution<float> diceroll(0.0f, 1.0f); 
for(int i-0; i<N; i++) 
{ 
    a[i] = diceroll(randomGenerator); 
    b[i] = diceroll(randomGenerator); 
} 

Теперь предположим, что необходимо суммировать a и b и сохранить результат в c, который в скалярном формате выглядит следующим образом:

for(int i=0; i<N; i++) 
{ 
    c[i] = a[i] + b[i]; 
} 

Какова будет версия в формате SSE2 с указанным выше кодом, учитывая, что входы: a и b, как определено выше (т. в качестве коллекции float), а ehe - c (также коллекция float)?

После довольно много исследований, я смог придумать следующее:

for(int i=0; i<N; i+=4) 
{ 
    float a_toload[4] = { a[i], a[i + 1], a[i + 2], a[i + 3] }; 
    float b_toload[4] = { b[i], b[i + 1], b[i + 2], b[i + 3] }; 
    __m128 loaded_a = _mm_loadu_ps(a_toload); 
    __m128 loaded_b = _mm_loadu_ps(b_toload); 

    float result[4] = { 0, 0, 0, 0 }; 
    _mm_storeu_ps(result, _mm_add_ps(loaded_a , loaded_b)); 
    c[i] = result[0]; 
    c[i + 1] = result[1]; 
    c[i + 2] = result[2]; 
    c[i + 3] = result[3]; 
} 

Однако, это, кажется, действительно громоздко и, конечно, весьма неэффективно: версия SIMD выше, фактически в три раза медленнее, чем исходная скалярная версия (измеренная, конечно, с оптимизацией в режиме выпуска Microsoft VS15 и после 1 миллиона итераций, а не только 12).

Answer 1

Ваш цикл может быть упрощено до

const int aligendN = N - N % 4; 
for (int i = 0; i < alignedN; i+=4) { 
    _mm_storeu_ps(&c[i], 
        _mm_add_ps(_mm_loadu_ps(&a[i]), 
        _mm_loadu_ps(&b[i]))); 
} 
for (int i = alignedN; i < N; ++i) { 
    c[i] = a[i] + b[i]; 
} 

Некоторые дополнительные пояснения:
1, небольшой обработки последних нескольких поплавков цикл бросить общий и когда N%4 != 0 или N неизвестно во время компиляции является обязательным ,
2, я замечаю, что вы выбираете нестандартную загрузку/сохранение версии, есть небольшой штраф по сравнению с выровненной версией. Я нашел эту ссылку в stackoverflow: Is the SSE unaligned load intrinsic any slower than the aligned load intrinsic on x64_64 Intel CPUs?

Answer 2

Вам не нужны промежуточные массивы для загрузки в регистры SSE. Просто загрузите непосредственно из своих массивов.

auto loaded_a = _mm_loadu_ps(&a[i]); 
auto loaded_b = _mm_loadu_ps(&b[i]); 
_mm_storeu_ps(&c[i], _mm_add_ps(loaded_a, loaded_b)); 

Вы также можете опустить два loaded переменные и их включения в оных, хотя компилятор должен сделать это для вас.

Вы должны быть осторожны с этим, так как это не сработает правильно, если размеры векторов не кратные 4 (вы получите доступ к концу массива, что приведет к неопределенному поведению и записи прошлого конец c может повредить).