Как векторизовать расчет расстояний с использованием SSE2

Question

A и B - векторы или длина N, где N может находиться в диапазоне от 20 до 200. Я хочу рассчитать квадрат расстояния между этими векторами, i.e. d^2 = || A-B ||^2.

До сих пор у меня есть:

float* a = ...; 
float* b = ...; 
float d2 = 0; 

for(int k = 0; k < N; ++k) 
{ 
    float d = a[k] - b[k]; 
    d2 += d * d; 
} 

Это, кажется, работает хорошо, за исключением того, что я профилированного мой код, и это узкое место (более 50% времени тратится только делает это). Я использую Visual Studio 2012, на Win 7, с этими вариантами оптимизации: /O2 /Oi /Ot /Oy-. Я понимаю, что VS2012 должен автоиндексировать этот цикл (используя SSE2). Однако, если я вставляю #pragma loop(no_vector) в код, я не получаю заметное замедление, поэтому я думаю, что цикл не является векторизованным. Компилятор подтверждает, что с этим сообщением:

info C5002: loop not vectorized due to reason '1105' 

Мои вопросы:

1. Можно ли исправить этот код, чтобы VS2012 может векторизации это?
2. Если нет, было бы целесообразно попытаться сами векторизовать код?
3. Можете ли вы порекомендовать мне веб-сайт, чтобы узнать о кодировании SSE2?
4. Есть ли какое-то значение N, ниже которого векторизация будет контрпродуктивной?
5. Что такое reason '1105'?

Answer 1

Из MSDN documentation код ошибки 1105 означает, что компилятор не может определить, как уменьшить код до векторизованных инструкций. Для операций с плавающей запятой указывается, что вам нужно указать параметр/fp: fast, чтобы включить любые сокращения с плавающей запятой.

Answer 2

Это довольно просто реализовать это с помощью SSE встроенных функций:

#include "pmmintrin.h" 

__m128 vd2 = _mm_set1_ps(0.0f); 
float d2 = 0.0f; 
int k; 

// process 4 elements per iteration 
for (k = 0; k < N - 3; k += 4) 
{ 
    __m128 va = _mm_loadu_ps(&a[k]); 
    __m128 vb = _mm_loadu_ps(&b[k]); 
    __m128 vd = _mm_sub_ps(va, vb); 
    vd = _mm_mul_ps(vd, vd); 
    vd2 = _mm_add_ps(vd2, vd); 
} 

// horizontal sum of 4 partial dot products 
vd2 = _mm_hadd_ps(vd2, vd2); 
vd2 = _mm_hadd_ps(vd2, vd2); 
_mm_store_ss(&d2, vd2); 

// clean up any remaining elements 
for (; k < N; ++k) 
{ 
    float d = a[k] - b[k]; 
    d2 += d * d; 
} 

Обратите внимание, что если вы можете гарантировать, что a и b являются 16 байт выровнены, то вы можете использовать _mm_load_ps, а не _mm_loadu_ps, которые могут помочь производительности, особенно на более старых (до Nehalem) процессорах.

Обратите внимание, что для таких циклов, где это очень мало арифметических инструкций относительно количества нагрузок, производительность может быть ограничена полосой пропускания памяти, и ожидаемое ускорение от векторизации может не реализоваться на практике.