Как получить 2 двойных или 4 числа с плавающей запятой, используя набор инструкций SSE? (До SSE4) - PullRequest
16 голосов
/ 01 апреля 2011

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

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 или, возможно, это переворачивание знака - ТАК однако слабость второго потребует условной проверки.

Ответы [ 3 ]

37 голосов
/ 13 мая 2011

Предлагаю использовать поразрядно и с маской.Положительные и отрицательные значения имеют одинаковое представление, отличается только старший значащий бит, это 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 теперь независимы), вы позволяете процессору выполнять сложения, которыепорядок.Поскольку инструкция конвейеризируется на современном процессоре, процессор может начать работу над новым дополнением до того, как закончится предыдущий.Кроме того, побитовые операции выполняются на отдельном исполнительном блоке, центральный процессор может фактически выполнять его в том же цикле, что и сложение / вычитание.Я предлагаю Руководства по оптимизации Agner Fog .

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

9 голосов
/ 13 мая 2011

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

x = _mm_max_ps(_mm_sub_ps(_mm_setzero_ps(), x), x)
2 голосов
/ 01 апреля 2011

Вероятно, самый простой способ заключается в следующем:

__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. Однако, если вы можете снизить точность до одинарной, то вы можете получить двукратное улучшение пропускной способности.

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

...