Можно ли векторизовать myNum + = a [b [i]] * c [i]; на x86_64?

Question

Какие свойства я бы использовал для векторизации следующего (если это вообще возможно) в x86_64?

double myNum = 0;
for(int i=0;i<n;i++){
    myNum += a[b[i]] * c[i]; //b[i] = int, a[b[i]] = double, c[i] = double
}

Answer 1

Вот мой пример, полностью оптимизирован и протестирован:

#include <emmintrin.h>

__m128d sum = _mm_setzero_pd();
for(int i=0; i<n; i+=2) {
    sum = _mm_add_pd(sum, _mm_mul_pd(
        _mm_loadu_pd(c + i),
        _mm_setr_pd(a[b[i]], a[b[i+1]])
    ));
}

if(n & 1) {
    sum = _mm_add_pd(sum, _mm_set_sd(a[b[n-1]] * c[n-1]));
}

double finalSum = _mm_cvtsd_f64(_mm_add_pd(
    sum, _mm_shuffle_pd(sum, sum, _MM_SHUFFLE2(0, 1))
));

Это дает очень красивый ассемблерный код с использованием gcc -O2 -msse2 (4.4.1).

Как вы можете сказать, наличие четного n сделает этот цикл более быстрым, чем выравнивание c. Если вы можете выровнять c, измените _mm_loadu_pd на _mm_load_pd для еще более быстрого выполнения.

Answer 2

Я бы начал с развертывания цикла. Что-то вроде

double myNum1 = 0, myNum2=0;
for(int i=0;i<n;i+=2)
{
    myNum1 += a[b[ i ]] * c[ i ];
    myNum2 += a[b[i+1]] * c[i+1];
}
// ...extra code to handle the remainder when n isn't a multiple of 2...
double myNum = myNum1 + myNum2;

Надеемся, что это позволяет компилятору чередовать нагрузки с арифметикой; профиль и посмотрите на сборку, чтобы увидеть, есть ли улучшения. В идеале компилятор будет генерировать инструкции SSE, но я не буду, если это произойдет на практике.

Дальнейшее развертывание может позволить вам сделать это:

__m128d sum0, sum1;
// ...initialize to zero...
for(int i=0;i<n;i+=4)
{
    double temp0 = a[b[ i ]] * c[ i ];
    double temp1 = a[b[i+1]] * c[i+1];
    double temp2 = a[b[i+2]] * c[i+2];
    double temp3 = a[b[i+3]] * c[i+3];
    __m128d pair0 = _mm_set_pd(temp0, temp1);
    __m128d pair1 = _mm_set_pd(temp2, temp3);
    sum0 = _mm_add_pd(sum0, pair0);
    sum1 = _mm_add_pd(sum1, pair1);
}
// ...extra code to handle the remainder when n isn't a multiple of 4...
// ...add sum0 and sum1, then add the result's components...

(извинения за псевдокод в начале и в конце, я думаю, что важной частью был цикл). Я не знаю точно, будет ли это быстрее; это зависит от различных задержек и от того, насколько хорошо компилятор может все переставить. Убедитесь, что вы профиль до и после, чтобы увидеть, было ли реальное улучшение.

Надеюсь, это поможет.

Answer 3

Процессоры Intel могут выполнять две операции с плавающей запятой, но одну загрузку за такт, поэтому доступ к памяти является самым жестким ограничением. Имея это в виду, я вначале стремился использовать упакованные нагрузки, чтобы уменьшить количество инструкций загрузки, и использовал упакованную арифметику только потому, что это было удобно. С тех пор я понял, что насыщение полосы пропускания памяти может быть самой большой проблемой, и весь беспорядок с инструкциями SSE мог бы быть преждевременной оптимизацией, если бы целью было ускорить выполнение кода, а не учиться векторизации.

SSE

Наименьшее количество возможных нагрузок без предположения о показателях в b требует развертывания цикла четыре раза. Одна 128-битная загрузка получает четыре индекса от b, две 128-битные загрузки каждый получает пару смежных двойных чисел от c, а сбор a требует независимых 64-битных нагрузок. Это пол 7 циклов на четыре итерации для последовательного кода. (достаточно для насыщения пропускной способности моей памяти, если доступ к a не кешируется). Я пропустил некоторые раздражающие вещи, такие как обработка нескольких итераций, которые не кратны 4.

entry: ; (rdi,rsi,rdx,rcx) are (n,a,b,c)
  xorpd xmm0, xmm0
  xor r8, r8
loop:
  movdqa xmm1, [rdx+4*r8]
  movapd xmm2, [rcx+8*r8]
  movapd xmm3, [rcx+8*r8+8]
  movd   r9,   xmm1
  movq   r10,  xmm1
  movsd  xmm4, [rsi+8*r9]
  shr    r10,  32
  movhpd xmm4, [rsi+8*r10]
  punpckhqdq xmm1, xmm1
  movd   r9,   xmm1
  movq   r10,  xmm1
  movsd  xmm5, [rsi+8*r9]
  shr    r10,  32
  movhpd xmm5, [rsi+8*r10]
  add    r8,   4
  cmp    r8,   rdi
  mulpd  xmm2, xmm4
  mulpd  xmm3, xmm5
  addpd  xmm0, xmm2
  addpd  xmm0, xmm3
  jl loop

Получение индексов - самая сложная часть. movdqa загружает 128 бит целочисленных данных из 16-байтового выровненного адреса (Nehalem имеет штрафы за задержку за смешивание инструкций SSE "integer" и "float"). punpckhqdq перемещает старшие 64 бита в младшие 64 бита, но в целочисленном режиме, в отличие от более простого имени movhlpd. 32-битные сдвиги выполняются в регистрах общего назначения. movhpd загружает один дубль в верхнюю часть регистра xmm, не нарушая нижнюю часть - это используется для загрузки элементов a непосредственно в упакованные регистры.

Этот код заметно быстрее, чем приведенный выше код, который, в свою очередь, быстрее, чем простой код, и для каждого шаблона доступа, но в простом случае B[i] = i, где наивный цикл на самом деле самый быстрый. Я также попробовал кое-что, например, функцию около SUM(A(B(:)),C(:)) в Фортране, которая в итоге оказалась эквивалентной простому циклу.

Я тестировал на Q6600 (65 нм Core 2 с частотой 2,4 ГГц) с 4 ГБ памяти DDR2-667 в 4 модулях. Тестирование пропускной способности памяти дает около 5333 МБ / с, поэтому кажется, что я вижу только один канал. Я собираю с Debian gcc 4.3.2-1.1, -O3 -Ffast-math -msse2 -Ftree-vectorize -std = gnu99.

Для тестирования я принимаю n равным одному миллиону, инициализируя массивы так, чтобы a[b[i]] и c[i] оба равнялись 1.0/(i+1), с несколькими различными шаблонами индексов. Один выделяет a с миллионом элементов и устанавливает b для случайной перестановки, другой выделяет a с 10M элементами и использует каждые 10-е, а последний выделяет a с 10M элементами и устанавливает b[i+1], добавляя случайное число от 1 до 9 b[i]. Я рассчитываю, сколько времени занимает один вызов с gettimeofday, очищая кеши, вызывая clflush над массивами, и измеряя 1000 испытаний каждой функции. Я построил сглаженные дистрибутивы времени исполнения, используя некоторый код из критерия (в частности, оценки плотности ядра в пакете statistics).

Пропускная способность

Теперь, для фактического важного примечания о пропускной способности. 5333 МБ / с с тактовой частотой 2,4 ГГц составляет чуть более двух байтов за цикл. Мои данные достаточно длинные, чтобы их нельзя было кешировать, и умножение времени выполнения моего цикла на (16 + 2 * 16 + 4 * 64) байт, загруженных за одну итерацию, если все пропускает, дает мне почти точно пропускную способность ~ 5333 МБ / с, которую имеет моя система , Должно быть довольно легко насытить эту пропускную способность без SSE. Даже если предположить, что a было полностью кэшировано, просто чтение b и c за одну итерацию перемещает 12 байтов данных, и наивный может начать новую итерацию в третьем цикле с конвейерной обработкой.

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

Более широкое оборудование может иметь большее значение. Система Nehalem, работающая с тремя каналами DDR3-1333, должна будет перемещать 3 * 10667 / 2,66 = 12,6 байта за такт для насыщения полосы пропускания памяти. Это было бы невозможно для одного потока, если a помещается в кеш - но при 64 байтах пропускается кэш строки в векторе, который быстро складывается - только одна из четырех нагрузок в моем цикле, отсутствующих в кешах, поднимает среднюю требуемую пропускную способность для 16 байт / цикл.

Answer 4

Это не будет векторизацией, как есть, из-за двойной косвенности индексов массива. Поскольку вы работаете с двойными кодами, SSE практически ничего не выиграет, тем более что большинство современных процессоров в любом случае имеют 2 FPU.

Answer 5

короткий ответ нет. Длинный ответ да, но не качественно. Вы будете подвергнуты штрафу за выполнение невыровненных нагрузок, что сведет на нет любую выгоду. Если вы не можете гарантировать, что b [i] последовательные индексы выровнены, вы, скорее всего, будете иметь худшую производительность после векторизации

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

чтобы ответить на ваш комментарий, вы должны сосредоточиться на массиве. Проще всего сразу попробовать заблокировать цикл в два раза, загрузить низко и высоко отдельно, а затем использовать мм * _ pd, как обычно. Псевдокод:

__m128d a, result;
for(i = 0; i < n; i +=2) {
  ((double*)(&a))[0] = A[B[i]];
  ((double*)(&a))[1] = A[B[i+1]];
  // you may also load B using packed integer instruction
  result = _mm_add_pd(result, _mm_mul_pd(a, (__m128d)(C[i])));
}

Я не помню точно имена функций, возможно, стоит перепроверить. Кроме того, используйте ключевое слово restrict с указателями, если вы знаете, что не может быть проблем с алиасами. Это позволит компилятору быть намного более агрессивным.