Сверхскоростное перемножение массивов

[Andrey_Kireew]

Как то на досуга, анализирую один программный код, с удивлением обнаружил (с помощью clock()) что основная его трудоёмкость (почти 98%) сосредоточена в одном месте - там, где выполняется перемножение векторов.
Операция выполняется такая $a[ii]-=b[ii]*w;$. Лично я считаю - проблема в том, что a и b одновременно к кэш никак не помещаются (их размеры примерно по 16 КБайт), и константа w тоже наверное какую то роль играет. В общем, думаю, что тормозит от многократных перезаписей в кэш.

У меня есть идея, перед выполнением операции, скопировать a и b в один буфер, поочерёдно, через 1 элемент. Константу w наверное нужно всегда держать вначале этого буфера. Так соседние элементы всегда будут рядом и многократных перегрузок кэш удастся избежать. Но тут возникают накладные расхода на загрузку в буфер, и выгрузку из него.

Может есть у кого ещё какие нибудь идеи, по поводу того как можно оптимизировать эту операцию?

[=SSN=]

Andrey_Kireew в сообщении #1443898 писал(а):

У меня есть идея, перед выполнением операции, скопировать a и b в один буфер, поочерёдно, через 1 элемент. Константу w наверное нужно всегда держать вначале этого буфера. Так соседние элементы всегда будут рядом и многократных перегрузок кэш удастся избежать.

Этим вы ничего не добьетесь. Оба массива по факту будут храниться в разных строках кэша. Количество элементов из каждого массива зависит от длины строки кэша.
Для ускорения перемножения матриц нужен процессор с поддержкой команд Intel® AVX-512 и компилятор, который понимает такие команды.

[mihaild]

Как раз AVX512 не факт что нужно. Но в целом убедитесь в том, что компилятор векторизует все операции, или векторизуйте их вручную, или воспользуйтесь библиотекой, где уже векторизованы. В blas нужная функция называется ?axpy.

[Dmitriy40]

Для ориентира: на Вашем CPU (intel core i7-6700) скорость вычислений может быть в пределе до 60 млрд итераций в секунду (для double, для single ещё вдвое быстрее), если использовать AVX2 команды FMA (4 итерации за такт на каждом ядре).
Но вероятно на GPU скорость будет ещё раз в 10-30 быстрее, тут не глядя сложно оценить (я не умею). Вы три года назад занимались перемножением матриц, тут очень похоже и даже проще.
Копировать в буфер смысла не имеет, процессор вполне поддерживает аппаратную предвыборку из двух массивов, а запись тормозить не должна.
Ну и оба массива точно помещаются в кэш L1 (у Вас он по 64КБ) каждого ядра.

[Andrey_Kireew]

Dmitriy40 в сообщении #1443960 писал(а):

скорость вычислений может быть в пределе до 60 млрд итераций в секунду (для double

по факту около 2 млрд. оп./сек, видимо это слишком мало

[Andrey_Kireew]

Уточните пожалуйста Dmitriy40, правильно ли я понимаю, что если вычисления идут на 1 ядре и без всяких там AVX2 команд FMA, а по старинке - то приведённую Вами цифру нужно делить на 4х4=16, т.е. получится что то около 3,75 млрд. оп./ сек, что довольно близко к наблюдаемому?

[mihaild]

Dmitriy40, а почему 4? У _mm256_fmadd_pd CPI 0.5.

[Dmitriy40]

mihaild
Действительно, спасибо, почему-то я помнил что блок FMA один, а не два, видимо из-за скалярных умножений MUL/IMUL/MULX, в основном ими пользуюсь, на ширину AVX домножил, а про второй блок не вспомнил, и в pdf поленился глянуть.
Но! В процессоре лишь один порт записи результатов в память, следовательно можно выполнить не более 4 double итераций за такт (сохранить лишь один регистр). И лишь два порта чтения памяти, потому прочитать можно не более двух регистров за такт, например a[] и b[], что тоже ограничивает до 4 double итераций на такт. Потому коэффициент таки 4 вместо 8. И выигрыша от применения FMA вместо AVX нет.
Далее все цифры уже строго по pdf на Skylake от Агнера Фога.

Andrey_Kireew в сообщении #1444026 писал(а):

если вычисления идут на 1 ядре и без всяких там AVX2 команд FMA, а по старинке - то приведённую Вами цифру нужно делить на 4х4=16, т.е. получится что то около 3,75 млрд. оп./ сек, что довольно близко к наблюдаемому?

Если по старинке, т.е. на FPU, то FMUL занимает 1 такт в потоке и выполняется в 0 порту, команда FSUB выполняется за 1 такт в потоке в 5 порту, две команды FLD выполняются параллельно в любом из 2 и 3 портов за 1 такт в потоке, плюс команда FSTP в 4+7 портах за 1 такт в потоке. Итого получается 1 операция обязательно в 0 порту, одна операция обязательно в 5 порту, две операции в любом из 2 или 3 портов и одна операция в 4+7 портах, в идеале должно уложиться в 1 такт во всех задействованных портах, т.е. 1 такт на итерацию. Надеюсь $w$ Вы/компилятор будете держать в FPU регистре. При частоте 3.4-4 ГГц это составит 3.4-4 млрд итераций в секунду (на ядро), что вдвое выше измеренного Вами.

С другой стороны, вместо FPU компиляторы давно перешли на использование xmm регистров даже для скалярных вычислений, посчитаем для них скорость: MULSD и SUBSD выполняются в любом из 0 или 1 портах и занимают 1 такт в потоке, MOVSD в обе стороны выполняются в других портах (2/3 и 4+7) и занимают тот же 1 такт в потоке. Итого все нужные операции могут выполниться за 1 такт (в портах 0 (MULSD), 1 (SUBSD), 2 (MOVSD из памяти), 3 (MULSD из памяти), 4+7 (MULSD в память), всего 5 операций за такт при максимум 6-ти допустимых) и дать скорость равную частоте процессора в 3.4-4 млрд итераций в секунду (на ядро).

Дополнительно оценим влияние команд организации цикла: DEC/SUB плюс JNZ спарятся и могут выполниться в свободном от вычислений 6-м порту за 1 такт, что на скорость не повлияет никак.

Раз у Вас скорость вдвое меньше, значит в поток вмешивается ещё какая-то лишняя в идеале команда, например загрузка в регистр величины $w$ или счётчика цикла. Или PREFETCH в каждом цикле. Или, скорее, и то и другое, т.к. по отдельности они займут лишь по полтакта в потоке и снизят скорость в полтора раза вместо двух. Хотя если счётчик цикла хранится в памяти (что идиотизм), то снижение скорости будет как раз ровно в два раза. Или Вы немного ошиблись с измерением скорости и она порядка 2.3 млрд итераций в секунду на частоте 3.4 ГГц (хотя у меня под даже ещё большей нагрузкой частота повышается на два шага, у Вас должна до 3.6-3.7 ГГц и 2.4 млрд итераций в секунду). На ядро разумеется.

Это всё было для скалярных вычислений, без векторизации цикла.
С векторизацией (руками или компилятором) скорость должна быть порядка 15 млрд итераций в секунду на ядро (частоту умножить на 1 операцию MOVDQA в память за такт и на 4 итерации в регистре для double). В памяти производится три операции (два чтения и одна запись) на 4 итерации, потому пропускная способность памяти тоже может ограничивать скорость, она должна быть не менее 24х скорости итераций (для double). Для распространённой памяти DDR4 2133МГц в двухканальном режиме полоса пропускания составляет 34 ГБ/с, что эквивалентно всего лишь 1.4 млрд double итераций в секунду (векторно!). Это применимо для объёмов данных примерно втрое и более больше объёма кэша L3, меньшие объёмы будут эффективно кэшироваться.

По идее указанные скорости не должны зависеть от размещения массивов в памяти или кэше данных любого уровня (за одним исключением для векторизированного цикла, см. выше) так как процессор аппаратно поддерживает предвыборку двух массивов из памяти и к моменту запроса данных (для массивов достаточной длины, более сотен элементов, чтобы латентность памяти незначительно повлияла на время вычислений) они уже окажутся в кэше L1. Для гарантии этого можно в код добавить команду PREFETCH, которая хоть и займёт любой из портов 2 или 3, но её можно подавать лишь дважды (для двух массивов) для 64 байтов линии кэша L1 (или 8 итераций), т.е. в идеале она замедлит выполнение лишь на $1\cdot2/8/2=13\%$ (такты∙команды/итераций/портов), что явно несущественно.

-- 10.03.2020, 20:10 --

Andrey_Kireew в сообщении #1444026 писал(а):

т.е. получится что то около 3,75 млрд. оп./ сек, что довольно близко к наблюдаемому?

Простите, не удержался, это не "довольно близко", это аж вдвое! "Довольно близко" это процентов 10-20, ну пусть даже 25, столько я обычно готов списать на точность измерений. Если у Вас точность измерений составляет разы, то ... я зря распинаюсь про детали вычислений в процессоре, уж "сравнимо в разы" будет выполняться почти любой код. Мне не жалко, но смысл теряется. Ещё раз простите.

[Andrey_Kireew]

Dmitriy40 в сообщении #1444089 писал(а):

тите, не удержался, это не "довольно близко", это аж вдвое! "Довольно близко" это процентов 10-20

В данном случае для меня 2 раза - вполне нормально, в том плане, что лучше сосредоточится на главной задаче, чем распыляться по мелочам. Потом, при желании, к этому всегда можно будет вернуться. Тем более, что алгоритм ещё сырой, впоследствии - может сильно измениться.

[feedinglight]

Опции компилятора пробовали выставлять?
Для gcc -march=native -O3
Для msvc вроде /arch:..., может ещё что-то.
Если компилятор достаточно не свежий, то он может и не знать про новые инструкции.

[Andrey_Kireew]

feedinglight в сообщении #1444613 писал(а):

Если компилятор достаточно не свежий, то он может и не знать про новые инструкции

vc 2010, наверное он такие инструкции не поддерживает, по крайней мере у меня ничего не вышло

[=SSN=]

Intel® C++ Compiler:

Цитата:

- Advanced optimization, Single Instruction Multiple Data (SIMD) vectorization, and loop and memory transformations
- Integration with Intel® Performance Libraries
- Using the latest OpenMP* parallel programming models

[Andrey_Kireew]

Это всё видимо касается последней версии, у меня же vc 2010

[Dmitriy40]

2000 итераций даже при Вашей скорости в два миллиарда в секунду занимают порядка одной микросекунды. Вы уверены что тормозит именно этот кусок? У меня подозрение, что тормозит не вычисление произведения, а память или ещё что-то.

Если очень надо, могу написать функции vecsubmuld(void * a, void * b, int n, double w), vecsubmuls(void * a, void * b, int n, single w) на ассемблере с AVX (и FMA) и скомпилировать в obj или dll файл, подключите к своему проекту. Расчётное время выполнения 0.13-1.0мкс (меньше если данные уже в L1, для single от 0.07мкс). Но боюсь что вызов её съест значительную часть выигрыша скорости.

[Andrey_Kireew]

Да, vc2010 не поддерживает AVX и FMA, попробую поставить компилятор посвежее, наверное начну с gcc4.7