Наборы инструкций SIMD Intel для 2D Matrix
Я разрабатываю высокопроизводительные алгоритмы на основе наборов инструкций Intel (AVX, FMA, ...). Мои алгоритмы (мои ядра) работают довольно хорошо, когда данные хранятся последовательно. Однако теперь я столкнулся с большой проблемой, и я не нашел обходного пути или решения для нее: см. 2D Matrix
int x, y; x = y = 4096;
float data[x*y]__attribute__((aligned(32)));
float buffer[y]__attribute__((aligned(32)));
/* simple test data */
for (i = 0; i < x; i++)
for (j = 0; j < y; j++)
data[y*i+j] = y*i+j; // 0,1,2,3...4095, | 4096,4097, ... 8191 |...
/* 1) Extract the columns out of matrix */
__m256i vindex; __m256 vec;
vindex = _mm256_set_epi32(7*y, 6*y, 5*y, 4*y, 3*y, 2*y, y, 0);
for(i = 0; i < x; i+=8)
{
vec = _mm256_i32gather_ps (&data[i*y], vindex, 4);
_mm256_store_ps (buffer[i], vec);
}
/* 2) Perform functions */
fft(buffer, x) ;
/*3) write back buffer into matrix*/
/* strided write??? ...*/
Я хочу найти очень эффективный способ сделать следующее:
Извлеките столбцы из матрицы: col1 = 0, 4096, 8192, ... col2 = 1, 4097, 8193, ... Я попробовал это с collect_ps, который действительно медленный.
Выполните мои высокоэффективные алгоритмы на извлеченных столбцах...
- Сохраните обратно столбцы в матрицу:
Есть ли какой-то особенный трюк для этого? Как вы можете читать и писать с ходу (например, 4096), используя наборы инструкций Intel?
Или есть какой-либо вариант манипуляции с памятью, чтобы получить столбцы из матрицы?
Спасибо!
3 ответа
[Для основных данных строки SIMD-доступ к строке быстрый, но к столбцу медленный]
Да, это характер x86-64 и подобных архитектур. Доступ к последовательным данным в памяти быстрый, но доступ к разрозненным данным (случайным или регулярным образом) медленный. Это является следствием наличия кэшей процессора.
Существует два основных подхода: копировать данные в новый порядок, который упрощает шаблоны доступа, или выполнять вычисления в порядке, обеспечивающем лучшие шаблоны доступа.
Нет, нет никаких правил большого пальца или золотых уловок, которые заставляют все просто работать. На самом деле, даже сравнение различных реализаций является сложным, потому что существует очень много сложных взаимодействий (от задержек кэша до чередования операций, к кешу и шаблонам доступа к памяти), что результаты сильно зависят от конкретного оборудования и набора данных под рукой.
Давайте рассмотрим типичный пример случая, умножение матрицы на матрицу. Допустим, мы умножаем две матрицы 5×5 (c = a × b), используя стандартный порядок данных в мажорной строке C:
c00 c01 c02 c03 c04 a00 a01 a02 a03 a04 b00 b01 b02 b03 b04
c05 c06 c07 c08 c09 a05 a06 a07 a08 a09 b05 b06 b07 b08 b09
c10 c11 c12 c13 c14 = a10 a11 a12 a13 a14 × b10 b11 b12 b13 b14
c15 c16 c17 c18 c19 a15 a16 a17 a18 a19 b15 b16 b17 b18 b19
c20 c21 c22 c23 c24 a20 a21 a22 a23 a24 b20 b21 b22 b23 b24
Если мы запишем результат как вертикальные векторные регистры SIMD с пятью компонентами, мы получим
c00 a00 b00 a01 b05 a02 b10 a03 b15 a04 b20
c01 a00 b01 a01 b06 a02 b11 a03 b16 a04 b21
c02 = a00 × b02 + a01 × b07 + a02 × b12 + a03 × b17 + a04 × b22
c03 a00 b03 a01 b08 a02 b13 a03 b18 a04 b23
c04 a00 b04 a01 b09 a02 b14 a03 b19 a04 b24
c05 a05 b00 a06 b05 a07 b10 a08 b15 a09 b20
c06 a05 b01 a06 b06 a07 b11 a08 b16 a09 b21
c07 = a05 × b02 + a06 × b07 + a07 × b12 + a08 × b17 + a09 × b22
c08 a05 b03 a06 b08 a07 b13 a08 b18 a09 b23
c09 a05 b04 a06 b09 a07 b14 a08 b19 a09 b24
и так далее. Другими словами, если c имеет тот же порядок, что и bмы можем использовать регистры SIMD с последовательным содержимым памяти для обоих c а также bи только собрать a, Кроме того, SIMD регистрирует a У всех компонентов одинаковое значение.
Обратите внимание, однако, что b регистры повторяются для всех пяти рядов c, Так что, может быть, лучше инициализировать c к нулю, затем сделать дополнения с продуктами, имеющими такие же b SIMD регистры:
c00 a00 b00 c05 a05 b00 c10 a10 b00 c15 a15 b00 c20 a20 b00
c01 a00 b01 c06 a05 b01 c11 a10 b01 c16 a15 b01 c21 a20 b01
c02 += a00 × b02, c07 += a05 × b02, c12 += a10 × b02, c17 += a15 × b02, c22 += a20 × b02
c03 a00 × b03 c08 a05 b03 c13 a10 b03 c18 a15 b03 c23 a20 b03
c04 a00 × b04 c09 a05 b04 c14 a10 b04 c19 a15 b04 c24 a20 b04
Если мы перенесли a сначала вектор SIMD регистрируется для a также будет получать значения из последовательных областей памяти. На самом деле, если a достаточно большой, линеаризуя шаблон доступа к памяти для a тоже дает достаточно большой прирост скорости, чтобы быстрее делать транспонированную копию (используя uint32_t для поплавков и uint64_t для двоих; т.е. не использовать SIMD или с плавающей запятой вообще для транспонирования, а просто копировать хранилище в порядке транспонирования).
Обратите внимание, что ситуация с порядком данных по главному столбцу, то есть по порядку транспонирования данных по сравнению с описанным выше, очень похожа. Здесь глубокая симметрия. Например, если c а также b имеют одинаковый порядок данных, и a в обратном порядке данных вы можете эффективно SIMD-векторизовать матричный продукт без необходимости копировать какие-либо данные. Различается только суммирование, поскольку это зависит от порядка данных, а умножение матриц не является коммутативным (a×b!= B ×a).
Очевидно, основная проблема заключается в том, что векторные регистры SIMD имеют фиксированный размер, поэтому вместо использования полной строки в качестве регистра, как в примере выше, вы можете использовать только частичные строки. (Если число столбцов в результате не кратно ширине регистра SIMD, у вас также есть этот частичный вектор, о котором нужно беспокоиться.)
SSE и AVX имеют относительно большое количество регистров (8, 16 или 32, в зависимости от набора используемых расширений) и, в зависимости от конкретного типа процессора, могут выполнять некоторые векторные операции одновременно или, по крайней мере, с меньшим количеством задержки, если чередуются несвязанные векторные операции. Таким образом, даже выбор ширины блока для одновременной работы и того, является ли этот блок подобным расширенному вектору или больше похожим на блочную подматрицу, зависит от обсуждения, тестирования и сравнения.
Итак, как мне сделать матрично-матричное умножение наиболее эффективно с использованием SIMD?
Как я уже сказал, это зависит от набора данных. Боюсь, нет простых ответов.
Основными параметрами (для выбора наиболее эффективного подхода) являются размеры и порядок в памяти мультипликатора и матрицы результатов.
Это становится еще интереснее, если вычислить произведение более двух матриц разных размеров. Это связано с тем, что количество операций зависит от порядка продуктов.
Почему вы так обескураживаете?
Я на самом деле нет. Все вышеперечисленное означает, что не так уж много людей могут справиться с подобными сложностями и оставаться в здравом уме и продуктивно, поэтому есть много неизученных подходов и много чего можно добиться в реальных условиях.
Даже если мы игнорируем встроенные компиляторы SIMD (<x86intrin.h> в этом случае), мы можем применить приведенную выше логику при проектировании внутренних структур данных, чтобы компилятор C, который мы используем, имел лучшие возможности для векторизации вычислений для нас. (Они пока не очень хороши в этом. Как я уже сказал, сложные вещи. Некоторым нравится Fortran лучше, чем C, потому что его выражения и правила упрощают компиляторам Fortran их оптимизацию и векторизацию.)
Если бы это было просто или легко, решения были бы хорошо известны. Но это не так, потому что это не так. Но это не значит, что это невозможно или вне нашей досягаемости; все, что это означает, - то, что достаточно умные разработчики еще не приложили достаточно усилий, чтобы разгадать это.
Если вы можете запустить свои алгоритмы на 8 (или 16 1) столбцах параллельно, одна обычная загрузка AVX может собрать 8 столбцов данных в один вектор. Затем другая загрузка может получить следующий ряд из всех этих столбцов.
Это имеет то преимущество, что вам никогда не нужно перетасовывать данные внутри вектора; все чисто вертикально, и у вас есть последовательные элементы каждого столбца в разных векторах.
Если бы это было сокращение, подобное суммированию столбца, вы бы получили 8 результатов параллельно. Если вы обновляете столбцы по ходу работы, то вы пишете векторы результатов сразу для 8 столбцов вместо вектора из 8 элементов одного столбца.
Сноска 1: 16 float столбцы = 64 байта = 1 полная строка кэша = два вектора AVX или один вектор AVX512. Чтение / запись полных строк кэша за раз намного лучше, чем сокращение одного столбца за раз, хотя обычно это хуже, чем доступ к последовательным строкам кэша. Особенно, если ваш шаг больше, чем страница 4k, предварительная выборка HW может не очень хорошо закрепиться на нем.
Очевидно, что для этого убедитесь, что ваши данные выровнены на 64, при этом длина строки тоже будет кратна 64 байтам. Заполните концы рядов, если вам нужно.
Делать только 1 вектор AVX (половину строки кэша) за раз было бы плохо, если первая строка будет удалена из L1d, прежде чем вы вернетесь назад, чтобы прочитать 2-й 32-байтовый вектор столбцов 8..15.
Другие предостережения:
Псевдонимы 4k могут быть проблемой: хранилище, а затем загрузка по адресам, кратным 4 кБ, сразу не обнаруживаются как неперекрывающиеся, поэтому загрузка блокируется хранилищем. Это может значительно снизить степень параллелизма, который может использовать процессор.
Шаги 4k также могут привести к пропущенным конфликтам в кеше, если вы касаетесь множества строк с псевдонимом одного и того же набора. Таким образом, обновление данных на месте может по-прежнему приводить к отсутствию кэша для хранилищ, поскольку строки могут быть удалены после загрузки и обработки до того, как хранилище будет готово к фиксации. Скорее всего, это будет проблемой, если шаг строки имеет большую степень 2. Если это окажется проблемой, возможно, выделите больше памяти в этом случае и добавьте в конце строки неиспользуемые элементы, поэтому формат хранения никогда не будет имеет большую мощность в 2 ряда.
Предварительная выборка соседней строки в кэше L2 (процессоры Intel) может попытаться заполнить пару каждой строки, к которой вы прикоснулись, если есть свободная пропускная способность. Это может привести к удалению полезных данных, особенно если вы близки к псевдонимам и / или емкости L2. Но если вы не выходите за эти пределы, это, вероятно, хорошо и поможет, когда вы перейдете к следующим 16 столбцам.
Данные должны храниться одна строка за другой в памяти. Поскольку C на самом деле не волнует, является ли это массив или матрица, вы можете получить доступ к элементам с
for(int i=0;i<columncount;i++)
data[i*LENGTH + desired_column];
теперь вы можете хранить данные или даже лучше адреса, чтобы передать их вашей рабочей функции. Если вы берете адреса, значения в матрице изменятся, поэтому вам не нужно записывать их обратно.