Можно ли векторизовать 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

короткий ответ №. Длинный ответ да, но неэффективно. Вы понесете штраф за выполнение несвязанных нагрузок, которые будут отрицать любую выгоду. Если вы не можете гарантировать, что 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]))); 
} 

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

Answer 3

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

Answer 4

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

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 5

Процессоры 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). 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 МБ/с, поэтому кажется, что я вижу только один канал. Я компилирую с gcc 4.3.2-1.1, -O3 -Ffast-math -msse2 -Ftree-vectorize -std = gnu99 в Debian.

Для тестирования я позволяю 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).

Bandwidth

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