Как абсолютное 2 двойных или 4 поплавка с использованием набора инструкций SSE? (До SSE4)

Здесь пример кода C, который я пытаюсь ускорить с помощью SSE, два массива имеют длину 3072 элемента с удвоением, могут опуститься до float, если мне не нужна точность удвоений.

double sum = 0.0;

for(k = 0; k < 3072; k++) {
    sum += fabs(sima[k] - simb[k]);
}

double fp = (1.0 - (sum / (255.0 * 1024.0 * 3.0)));

В любом случае моя текущая проблема заключается в том, как сделать шаг fabs в регистре SSE для удвоения или поплавка, чтобы я мог сохранить весь расчет в регистрах SSE, чтобы он оставался быстрым, и я могу распараллелить все этапы путем частичного разворота этот цикл.

Здесь некоторые ресурсы, которые я нашел fabs() asm или, возможно, этот переворачивающий знак - SO, однако слабость второго требует условной проверки.

Ответ 1

Вероятно, самый простой способ:

__m128d vsum = _mm_set1_pd(0.0);        // init partial sums
for (k = 0; k < 3072; k += 2)
{
    __m128d va = _mm_load_pd(&sima[k]); // load 2 doubles from sima, simb
    __m128d vb = _mm_load_pd(&simb[k]);
    __m128d vdiff = _mm_sub_pd(va, vb); // calc diff = sima - simb
    __m128d vnegdiff = mm_sub_pd(_mm_set1_pd(0.0), vdiff); // calc neg diff = 0.0 - diff
    __m128d vabsdiff = _mm_max_pd(vdiff, vnegdiff);        // calc abs diff = max(diff, - diff)
    vsum = _mm_add_pd(vsum, vabsdiff);  // accumulate two partial sums
}

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

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

Ответ 2

Я предлагаю использовать побитовое и с маской. Положительные и отрицательные значения имеют одно и то же представление, только самый старший бит отличается, он равен 0 для положительных значений и 1 для отрицательных значений, см. формат двойной точности чисел. Вы можете использовать один из них:

inline __m128 abs_ps(__m128 x) {
    static const __m128 sign_mask = _mm_set1_ps(-0.f); // -0.f = 1 << 31
    return _mm_andnot_ps(sign_mask, x);
}

inline __m128d abs_pd(__m128d x) {
    static const __m128d sign_mask = _mm_set1_pd(-0.); // -0. = 1 << 63
    return _mm_andnot_pd(sign_mask, x); // !sign_mask & x
}

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

double norm(const double* sima, const double* simb) {
__m128d* sima_pd = (__m128d*) sima;
__m128d* simb_pd = (__m128d*) simb;

__m128d sum1 = _mm_setzero_pd();
__m128d sum2 = _mm_setzero_pd();
for(int k = 0; k < 3072/2; k+=2) {
    sum1 += abs_pd(_mm_sub_pd(sima_pd[k], simb_pd[k]));
    sum2 += abs_pd(_mm_sub_pd(sima_pd[k+1], simb_pd[k+1]));
}

__m128d sum = _mm_add_pd(sum1, sum2);
__m128d hsum = _mm_hadd_pd(sum, sum);
return *(double*)&hsum;
}

Развертывая и разбивая зависимость (sum1 и sum2 теперь независимы), вы позволяете процессору выполнять дополнения нашего порядка. Поскольку инструкция конвейерна на современном процессоре, процессор может начать работу над новым дополнением до завершения предыдущего. Кроме того, побитовые операции выполняются на отдельном исполнительном устройстве, CPU может фактически выполнять его в том же цикле, что и сложение/вычитание. Я предлагаю руководства по оптимизации Agner Fog.

Наконец, я не рекомендую использовать openMP. Цикл слишком мал, и накладные расходы на распределение работы между несколькими потоками могут быть больше, чем любая потенциальная выгода.

Ответ 3

Максимум -x и x должен быть abs (x). Здесь он находится в коде:

x = _mm_max_ps(_mm_sub_ps(_mm_setzero_ps(), x), x)