Какие внутренние функции я бы использовал для векторизации следующего (если это вообще возможно) на 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
}
Какие внутренние функции я бы использовал для векторизации следующего (если это вообще возможно) на 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
}
Вот мой вариант, полностью оптимизированный и протестированный:
#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
для еще более быстрого выполнения.
finalSum = 0;
или finalsum = a[b[ n-1 ]] * c[ n-1 ];
, параллельно с горизонтальной суммой.
- person Peter Cordes; 06.10.2020
Я бы начал с разматывания петли. Что-то типа
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...
(извините за псевдокод в начале и в конце, я полагаю, что важной частью был цикл). Я не знаю точно, будет ли это быстрее; это зависит от различных задержек и от того, насколько хорошо компилятор может все переставить. Убедитесь, что вы профилируете до и после, чтобы увидеть, было ли реальное улучшение.
Надеюсь, это поможет.
Процессоры Intel могут выполнять две операции с плавающей запятой, но одну загрузку за цикл, поэтому доступ к памяти является самым жестким ограничением. Имея это в виду, я стремился сначала использовать упакованные загрузки, чтобы уменьшить количество инструкций загрузки, и использовал упакованную арифметику просто потому, что это было удобно. С тех пор я понял, что насыщение пропускной способности памяти может быть самой большой проблемой, и все возни с инструкциями 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). 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 байт. /цикл.
movd/q reg, xmm
. Более поздние процессоры (начиная с Haswell в ~ 2013 году) имеют для этого только пропускную способность 1 / такт, поэтому это будет узким местом. agner.org/optimize/instruction_tables.pdf (так что вам нужно сделать 64 -битовая скалярная загрузка для распаковки с помощью mov
/shr
). Но да, это наверное хорошо для Core2/Nehalem.
- person Peter Cordes; 06.10.2020
[rcx+8*r8+8]
- это должно быть +16
, иначе оно сместится и перекроет предыдущий movapd
.
- person Peter Cordes; 06.10.2020
короткий ответ нет. Длинный ответ да, но неэффективно. Вы понесете штраф за невыровненные нагрузки, которые сведут на нет все преимущества. Если вы не можете гарантировать, что последовательные индексы b[i] выровнены, у вас, скорее всего, будет худшая производительность после векторизации.
если вы заранее знаете, что такое индексы, лучше всего развернуть и указать явные индексы. Я сделал что-то подобное, используя специализацию шаблонов и генерацию кода. если интересно могу поделиться
чтобы ответить на ваш комментарий, вам в основном нужно сосредоточиться на массиве. Проще всего сразу попробовать заблокировать цикл в два раза, загрузить низкий и высокий уровни a отдельно, а затем использовать mm*_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])));
}
Я точно не помню названия функций, может стоит перепроверить. Кроме того, используйте ключевое слово ограничения с указателями, если вы знаете, что не может быть проблем с псевдонимами. Это позволит компилятору быть более агрессивным.
Это не будет векторизоваться как есть из-за двойной косвенности индексов массива. Поскольку вы работаете с двойниками, от SSE мало что можно получить, особенно потому, что большинство современных процессоров в любом случае имеют 2 FPU.
mulsd
и addsd
, почему бы не использовать одну и ту же инструкцию (которая стоит одинаковое количество циклов) для выполнения ДВОЙНОЙ работы? См. мой ответ о том, как наилучшим образом использовать SSE(2) в этом случае.
- person LiraNuna; 02.03.2010
_mm_loadu_pd
единственное настоящее узкое место, остальное почти ничто по сравнению с ним. Я рекомендовал выровнять c
и попытаться сохранить n
даже, но если это невозможно, то вы правы, этот код может работать с минимальной выгодой. Это все еще быстрее, чем оригинал, который использует слишком много addpd
и mulpd
.
- person LiraNuna; 02.03.2010