Это хорошая идея, но компиляторы уже сделают это для вас, если они знают, что это безопасно. Используйте double *restrict output и const double *restrict input, чтобы обещать компилятору, который хранит output[], не делатьизмените то, что будет считываться из input[].
Но автоматическая векторизация с SIMD является еще более важной оптимизацией , производящей 2 или 4 double результатов на инструкцию.GCC и ICC уже сделают это в -O3 после проверки на совпадение.(Но clang не может автоматически векторизовать это, просто развернув его с помощью скаляра [v]addsd, избегая ненужных перезагрузок.
К сожалению, ваша оптимизированная версия побеждает авто-векторизацию! (Это ошибка компилятора, т.е.ошибка пропущенной оптимизации, когда она знает, что вывод не перекрывается, поэтому перечитывание источника из памяти эквивалентно).
Похоже, что gcc довольно хорошо справляется с оригинальной версией, с -O3 -march=native (особенно при настройке на Intel, где более широкие векторы с AVX того стоят.) Я вычисляю 4 double результатов параллельно с 3-мя не выровненными нагрузками и 2 vaddpd ymm.
. Проверяетперекрытие перед использованием векторизованного цикла. Вы можете использовать double *restrict output и input, чтобы пообещать, что указатели не перекрываются, поэтому не требуется резервный цикл.
L1d кешПропускная способность на современных процессорах превосходна, перезагрузка одних и тех же данных не является большой проблемой (2 загрузки в такт) . Пропускная способность инструкций является большей проблемой. Источник памяти addsdне стоит намного дороже, чем хранение данных в регистрах.
Если векторизация с использованием 128-битных векторов имеет смысл сохранить вектор in[idx+1..2] для использования в качестве вектора in[idx+ -1..1] на следующей итерации.GCC фактически делает это.
Но когда вы создаете 4 результата на инструкцию, ни один из 3 входных векторов из одной итерации напрямую не пригодится для следующей.Хотя, возможно, было бы полезно сохранить некоторую полосу пропускания порта загрузки с помощью случайного перемешивания для создания одного из 3 векторов из результата загрузки.Я бы попробовал это, если бы я вручную векторизовал с __m256d встроенными функциями.Или с float с 128-битными __m128 векторами.
#define SIZE 1000000
void stencil_restrict(double *restrict input, double *restrict output)
{
int idx = 1;
output[0] = input[0] + input[1];
for(; idx < SIZE - 1; idx++){
output[idx] = input[idx-1] + input[idx] + input[idx+1];
}
output[idx] = input[idx-1] + input[idx];
}
компилируется в этот asm с gcc8.3 -O3 -Wall -std=c99 -march=broadwell -masm=intel, из проводника компилятора Godbolt (-ffast-math в этом случае не требуется и не имеет значения для внутреннего цикла.)
stencil_restrict:
vmovsd xmm0, QWORD PTR [rdi]
vaddsd xmm0, xmm0, QWORD PTR [rdi+8]
xor eax, eax
vmovsd QWORD PTR [rsi], xmm0 # first iteration
### Main loop
.L12:
vmovupd ymm2, YMMWORD PTR [rdi+8+rax] # idx +0 .. +3
vaddpd ymm0, ymm2, YMMWORD PTR [rdi+rax] # idx -1 .. +2
vaddpd ymm0, ymm0, YMMWORD PTR [rdi+16+rax] # idx +1 .. +4
vmovupd YMMWORD PTR [rsi+8+rax], ymm0 # store idx +0 .. +3
add rax, 32 # byte offset += 32
cmp rax, 7999968
jne .L12
# cleanup of last few elements
vmovsd xmm1, QWORD PTR [rdi+7999976]
vaddsd xmm0, xmm1, QWORD PTR [rdi+7999968]
vaddsd xmm1, xmm1, QWORD PTR [rdi+7999984]
vunpcklpd xmm0, xmm0, xmm1
vaddpd xmm0, xmm0, XMMWORD PTR [rdi+7999984]
vmovups XMMWORD PTR [rsi+7999976], xmm0
vmovsd xmm0, QWORD PTR [rdi+7999984]
vaddsd xmm0, xmm0, QWORD PTR [rdi+7999992]
vmovsd QWORD PTR [rsi+7999992], xmm0
vzeroupper
ret
К сожалению, gcc использует режимы индексированной адресации, поэтому инструкции vaddpd систочник памяти раскладывается в 2 мопа для внешнего интерфейса семейства SnB (включая ваш Broadwell Xeon E5-2698 v4). Режимы Micro Fusion и Addressing
vmovupd ymm2, YMMWORD PTR [rdi+8+rax] # 1 uop, no micro-fusion
vaddpd ymm0, ymm2, YMMWORD PTR [rdi+rax] # 2 uops. (micro-fused in decoders/uop cache, unlaminates)
vaddpd ymm0, ymm0, YMMWORD PTR [rdi+16+rax] # 2 uops. (ditto)
vmovupd YMMWORD PTR [rsi+8+rax], ymm0 # 1 uop (stays micro-fused, but can't use the port 7 store AGU)
add rax, 32 # 1 uop
cmp rax, 7999968 # 0 uops, macro-fuses with JNE
jne .L12 # 1 uop
Анализ пропускной способности, см. https://agner.org/optimize/ и Какие соображения относятся к прогнозированию задержки для операций на современных суперскалярных процессорах икак я могу вычислить их вручную?
Цикл GCC - это 8 мопов в слитых доменах для этапа выдачи / переименования переднего плана для отправки в неработающий бэкэнд.Это означает, что максимальная пропускная способность внешнего интерфейса составляет 1 итерацию на 2 цикла.
[v]addpd в Intel до того, как Skylake сможет работать только на порту 1, в отличие от [v]mulpd или FMA с удвоенной пропускной способностью.(Skylake отбросил выделенный модуль добавления FP и запускает добавление FP одинаково для mul и fma.) Так что это также узкое место с 2 циклами на итерацию.
У нас есть 3 загрузки + 1 хранилище, каждый из которых нуждается в одном изпорт 2 или 3. (Хранилища в режиме индексированной адресации не могут использовать AGU хранилища на порте 7).Так что это еще один 2 цикла на узкое место итерации.Но не совсем; невыровненные загрузки, которые пересекают границы строк кэша, стоят дороже.Эксперименты показывают, что Intel Skylake (и, вероятно, также Broadwell) воспроизводит загруженные мопы, которые, как обнаружено, разделены на строки кэша, поэтому они снова запускаются для получения данных из второй строки кэша. Как можно точно оценить скорость невыровненного доступа на x86_64 .
Наши данные выровнены по 8 байтов, но мы делаем 32-байтовые нагрузки, равномерно распределенные по всем 8-байтовым смещениям в пределах 64-байтовой строки.В 5 из этих 8 начальных элементов разделение строки кэша отсутствует.На остальных 3 есть.Таким образом, средняя стоимость действительно составляет 3 * (8+3)/8 = 4.125 загрузок мопов за одну итерацию.Я не знаю, нужно ли переиграть адреса магазина.Возможно нет;это важно, когда данные фиксируются из буфера хранилища в L1d, что имеет значение, а не для адреса хранилища или данных хранилища.(До тех пор, пока он не разделен по границе 4К, что будет произойдет с неправильно выровненным выводом).
Предполагается, что выходное выравнивание любого, кроме output[1], выровнено по 32 байта.Asm хранит output[0] вне цикла, а затем фактически делает output[i*4 + 1], поэтому каждое другое хранилище будет разделено строкой кэша.
В этом случае было бы лучше достичь границы выравнивания длявыходной массив.Gcc7 и более ранние версии предпочитают выравнивать один из указателей с помощью пролога цикла, но, к сожалению, в любом случае они выбирают вход, куда мы загружаем, из всех выравниваний.
В любом случае фактическим узким местом GCC является пропускная способность порта 2 / порта 3. В среднем 5,125 мопов за итерацию для этих 2 портов = теоретическая максимальная средняя пропускная способность 1 итерации (4 двойных) за 2,5625 цикла .
Использование неиндексированного хранилища будет иметьуменьшило это узкое место.
Но при этом игнорируются штрафы за разделение 4k, что составляет ~ 100 циклов на Бродвелле, и предполагается идеальная предварительная выборка HW, которая может поддерживать ~ 12,5 байт / цикл в каждом направлении (загружается и сохраняется). Таким образом, более вероятно, что это будет узким местом в пропускной способности памяти, если данные не были уже горячими в кэше L2. L1d может поглощать избыточные нагрузки тех же байтов, но все еще существует значительная не избыточная пропускная способность.
Небольшая развертка позволила бы незапланированному исполнению видеть впереди себя и помогать поглощать пузыри из-за промахов кэша, когда предварительная выборка HW не идет в ногу.Если бы он использовал неиндексированный режим адресации для хранилища, он мог бы использовать порт 7, уменьшая нагрузку на порты 2/3.Это позволило бы нагрузкам опережать добавления, надеясь поглотить пузырьки при пересечении
Повторное использование данных в регистрах с 128-битными векторами
Внутренний цикл из gcc8.3 -O3 -Wall -std=c99 -march=broadwell -mno-avx
# prologue to reach an alignment boundary somewhere?
.L12:
movupd xmm2, XMMWORD PTR [rdi+rax]
movupd xmm1, XMMWORD PTR [rdi+8+rax]
addpd xmm0, xmm2
addpd xmm0, xmm1
movups XMMWORD PTR [rsi+rax], xmm0
add rax, 16
movapd xmm0, xmm1 # x = z
cmp rax, 7999992
jne .L12
Это регрессия против gcc7.4, которая позволяет избежать копирования регистра.(Но gcc7 теряет накладные расходы цикла на счетчике, отдельном от индекса массива.)
# prologue to reach an alignment boundary so one load can be aligned.
# r10=input and r9=input+8 or something like that
# r8=output
.L18: # do {
movupd xmm0, XMMWORD PTR [r10+rdx]
add ecx, 1
addpd xmm0, xmm1 # x+y
movapd xmm1, XMMWORD PTR [r9+rdx] # z for this iteration, x for next
addpd xmm0, xmm1 # (x+y) + z
movups XMMWORD PTR [r8+rdx], xmm0
add rdx, 16
cmp ecx, r11d
jb .L18 # } while(i < max);
Это, вероятно, все же медленнее, чем 256-битные векторы AVX, в среднем.
С AVX для 128-битные векторы (например, настройка на Piledriver), он мог бы избежать отдельной нагрузки movupd xmm0 и использовать vaddpd xmm0, xmm1, [r10+rdx].
Оба они не используют выровненные хранилища, но также не могут воспользоватьсяскладывание загрузки в операнд памяти на addpd после нахождения известного выравнивания в input: /
Фактические эксперименты с производительностью на Skylake показывают, что реальная производительность довольно близка к той, которую я предсказывал, если данныепомещается в кэш L1d.
Интересный факт: со статическими буферами, такими как глобальные double in[SIZE+10];, gcc создаст версию цикла, которая использует неиндексированные режимы адресации.Это дает ускорение от ~ 800 мс до ~ 700 мс для многократного запуска его в цикле с SIZE = 1000.Более подробная информация будет обновлена позже.