Какие внутренние свойства я использовал бы для векторизации следующего (если это возможно, векторизации) на 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
на еще более быстрое время выполнения.
Я бы начал с разворачивания цикла. Что-то вроде
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 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(:))
в Fortran, которая в итоге была эквивалентна простому циклу.
Я тестировал на 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 проб каждой функции. Я построил сглаженные распределения времени выполнения, используя некоторый код из кишок criterion (в частности, оценка плотности ядра в пакете 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 байт/цикл.
короткий ответ №. Длинный ответ да, но неэффективно. Вы понесете штраф за выполнение несвязанных нагрузок, которые будут отрицать любую выгоду. Если вы не можете гарантировать, что b [i] последовательные индексы выровнены, у вас, скорее всего, будет хуже производительность после векторизации
если вы заранее знаете, какие индексы, ваше лучшее, что нужно развернуть и указать явные индексы. Я сделал что-то подобное, используя специализированную специализацию и генерацию кода. если вы заинтересованы, я могу поделиться
чтобы ответить на ваш комментарий, вам в основном нужно сосредоточиться на массиве. Самое простое, что можно попробовать сразу, это заблокировать цикл в два раза, загрузите минимум и максимум отдельно, а затем используйте 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.