Наборы инструкций SIMD Intel для 2D Matrix

Question

Я разрабатываю высокопроизводительные алгоритмы на основе наборов инструкций 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??? ...*/

Я хочу найти очень эффективный способ сделать следующее:

1. Извлечение столбцов из матрицы: col1 = 0, 4096, 8192, ... col2 = 1, 4097, 8193, ... Я попробовал это с collect_ps, который очень медленный.

2. Выполнить мои высокоэффективные алгоритмы на извлеченных столбцах ...

3. Сохраните обратно столбцы в матрицу:

Есть ли какой-нибудь особый трюк для этого? Как вы можете читать и писать с ходу (например, 4096), используя наборы инструкций Intel?

Или есть какая-нибудь опция манипулирования памятью для получения столбцов из матрицы?

Спасибо!

Answer 1

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

SSE и AVX имеют относительно большое количество регистров (8,16 или 32, в зависимости от набора используемых расширений) и в зависимости от конкретного типа процессора, может быть в состоянии выполнять некоторые векторные операции одновременно или, по меньшей мере, с меньшим временем ожидания, если чередуются несвязанные векторные операции.Таким образом, даже выбор ширины блока для одновременной работы и того, является ли этот блок похожим на расширенный вектор или более похожим на блочную подматрицу, зависит от обсуждения, тестирования и сравнения.

Так как мне сделать матрично-матричное умножение наиболее эффективно с использованием SIMD?

Как я уже сказал, это зависит от набора данных.Боюсь, нелегких ответов.

Основными параметрами (для выбора наиболее эффективного подхода) являются размеры и порядок в памяти матриц мультипликатора и результата.

Будет еще интереснее, если вы рассчитаете произведение более двух матриц разных размеров.Это потому, что количество операций зависит от порядка продуктов.

---

Почему вы так обескураживаете?

На самом деле это не так.Все вышеперечисленное означает, что не так уж много людей могут справиться с подобными сложностями и оставаться в здравом уме и продуктивно, поэтому существует множество неизученных подходов и много чего можно добиться в реальных условиях.

Даже еслимы игнорируем встроенные компиляторы SIMD (в данном случае <x86intrin.h>), мы можем применять приведенную выше логику при проектировании внутренних структур данных, так что используемый нами компилятор C имеет наилучшие возможности для векторизации вычислений для нас.(Они пока не очень хороши в этом. Как я уже сказал, сложные вещи. Некоторым нравится Fortran лучше, чем C, потому что его выражения и правила упрощают их компилятору и оптимизируют их.)

Если бы это было просто или легко, решения были бы хорошо известны.Но это не так, потому что это не так.Но это не значит, что это невозможно или вне нашей досягаемости;все, что это означает, - то, что достаточно умные разработчики еще не приложили достаточно усилий, чтобы разгадать это.

Answer 2

Если вы можете запустить свои алгоритмы для 8 (или 16 ¹ ) столбцов параллельно , одна обычная загрузка 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 столбцам.

Answer 3

Данные должны храниться одна строка за другой в памяти. Поскольку C на самом деле не волнует, является ли он массивом или матрицей, вы можете получить доступ к элементам с помощью

for(int i=0;i<columncount;i++) data[i*LENGTH + desired_column];

теперь вы можете хранить данные или, что еще лучше, адреса, чтобы передать их в вашу рабочую функцию. Если вы возьмете адреса, значения в матрице изменятся, поэтому вам не нужно записывать их обратно.