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

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


gcc sse
person Pharaun    schedule 01.04.2011    source источник

Ответы (3)


arrow_upward
2
arrow_downward

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

__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 в скалярное значение после цикла, но это довольно тривиально и не критично для производительности.

person Paul R    schedule 01.04.2011
comment
Во второй части, где я позаботился о частичных суммах, мне потребовалось немного времени, чтобы собраться вместе, но в целом решение сработало. Спасибо! - person Pharaun; 02.04.2011
comment
@Pharaun: спасибо за отзыв - вы сравнивали код SSE с исходным скалярным кодом, и если да, то насколько улучшилась производительность? - person Paul R; 02.04.2011
comment
@Paul с однопоточной реализацией я скинул 20%, но с многопоточной через openMP это был всего лишь процент, поэтому мне все еще нужно поработать над ее настройкой и изучить возможность перехода к числам с плавающей запятой для большей производительности. - person Pharaun; 03.04.2011
comment
@Paul, только что перевел его на реализацию с плавающей запятой, и это сократило время выполнения примерно на 15–30% по сравнению с простым C, и это даже с включенным openMP, так что, похоже, это помогло моей пропускной способности памяти + четырехъядерная перекачка числа с плавающей запятой. к крупной победе. - person Pharaun; 03.04.2011
comment
@Pharaun: отлично - рад, что это помогло. Если вы можете перейти на последний процессор Sandy Bridge и использовать AVX, вы должны увидеть гораздо большее улучшение. (AVX выполняет 8 поплавков по сравнению с 4 поплавками с SSE). - person Paul R; 03.04.2011
comment
@Paul Когда-нибудь, когда-нибудь! Но да, обновление до Sandy Bridge — это не то, что я могу сделать в данный момент. - person Pharaun; 03.04.2011
comment
И еще один анонимный голос против, 4 года спустя - неужели так трудно оставить комментарий, объясняющий, в чем проблема? Я рад улучшить или даже удалить ответ, если это необходимо, но я действительно не понимаю, что с ним не так? - person Paul R; 14.06.2015
comment
Это решение намного сложнее, чем должно быть, в частности, для его выполнения не нужны никакие арифметические действия. Переворот начального бита двоичного представления IEE754, как показано Норбертом П. ниже, — это путь. - person Wenzel Jakob; 03.08.2016
comment
@WenzelJakob: гораздо более сложный, кажется, значительно преувеличивает - вам все равно нужно вычесть два входных значения, поэтому единственная реальная разница в том, что я делаю вычитание и максимум в своем коде, а не немного переворачиваю Решение Норберта. Также возникает вопрос, нужно ли правильное поведение для крайних случаев, таких как NaNs — я не уверен, что решение с переворотом битов дает строго корректные результаты в таких случаях, но это может быть и не важно. В любом случае, точка зрения принята - спасибо, что нашли время, чтобы объяснить свое отрицательное голосование. - person Paul R; 04.08.2016

arrow_upward
37
arrow_downward

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

person Norbert P.    schedule 13.05.2011
comment
Когда у меня скоро будет время, я посмотрю на это, но я хотел прокомментировать аспект OpenMP... Фрагмент, который я использовал/представил в этом вопросе, был вложенным циклом внутри другого гораздо больший цикл, и я использовал OpenMP во внешнем цикле, а не во внутреннем цикле :) Но да, вы были бы правы, если бы OpenMP был для этого меньшего внутреннего цикла. - person Pharaun; 16.05.2011
comment
Привет, я смог реализовать это, и это ускорило работу на несколько процентов, что было приятно :) И мне очень понравились руководства Agner Fog, они были действительно хороши! - person Pharaun; 24.06.2011

arrow_upward
10
arrow_downward

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

x = _mm_max_ps(_mm_sub_ps(_mm_setzero_ps(), x), x)
person Christian Buchner    schedule 13.05.2011
comment
Да, я не знал об этом трюке, когда задавал вопрос, но да, это имеет смысл :) - person Pharaun; 16.05.2011