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

float compute_distance(const float* p1, const float* p2) const
{
__m256 euclidean = _mm256_setzero_ps();

const uint16_t n = nbr_dimensions;
const uint16_t aligend_n = n - n % 16;
const float* local_selected = selected_dimensions;

for (uint16_t i = 0; i < aligend_n; i += 16)
{
const __m256 r1 = _mm256_sub_ps(_mm256_load_ps(&p1[i]), _mm256_load_ps(&p2[i]));
euclidean = _mm256_fmadd_ps(_mm256_mul_ps(r1, r1), _mm256_load_ps(&local_selected[i]), euclidean);
const __m256 r2 = _mm256_sub_ps(_mm256_load_ps(&p1[i + 8]), _mm256_load_ps(&p2[i + 8]));
euclidean = _mm256_fmadd_ps(_mm256_mul_ps(r2, r2), _mm256_load_ps(&local_selected[i + 8]), euclidean);
}
float distance = hsum256_ps_avx(euclidean);

for (uint16_t i = aligend_n; i < n; ++i)
{
const float num = p1[i] - p2[i];
distance += num * num * local_selected[i];
}

return distance;
}

float compute_distance(const float* p1, const float* p2) const
{
const size_t n = nbr_dimensions;
const size_t aligend_n = n - n % 16;
const unsigned int* local_selected = selected_dimensions;
const __m256* local_masks = masks;

__m256 euc1 = _mm256_setzero_ps(), euc2 = _mm256_setzero_ps(),
euc3 = _mm256_setzero_ps(), euc4 = _mm256_setzero_ps();

const size_t n_max = aligend_n/8;
for (size_t i = 0; i < n_max; i += 4)
{
const __m256 r1 = _mm256_sub_ps(_mm256_load_ps(&p1[i * 8 + 0]), _mm256_load_ps(&p2[i * 8 + 0]));
const __m256 r1_1 = _mm256_and_ps(r1, local_masks[i + 0]);
euc1 = _mm256_fmadd_ps(r1_1, r1_1, euc1);

const __m256 r2 = _mm256_sub_ps(_mm256_load_ps(&p1[i * 8 + 8]), _mm256_load_ps(&p2[i * 8 + 8]));
const __m256 r2_1 = _mm256_and_ps(r2, local_masks[i + 1]);
euc2 = _mm256_fmadd_ps(r2_1, r2_1, euc2);

const __m256 r3 = _mm256_sub_ps(_mm256_load_ps(&p1[i * 8 + 16]), _mm256_load_ps(&p2[i * 8 + 16]));
const __m256 r3_1 = _mm256_and_ps(r3, local_masks[i + 2]);
euc3 = _mm256_fmadd_ps(r3_1, r3_1, euc3);

const __m256 r4 = _mm256_sub_ps(_mm256_load_ps(&p1[i * 8 + 24]), _mm256_load_ps(&p2[i * 8 + 24]));
const __m256 r4_1 = _mm256_and_ps(r4, local_masks[i + 3]);
euc4 = _mm256_fmadd_ps(r4_1, r4_1, euc4);
}

float distance = hsum256_ps_avx(_mm256_add_ps(_mm256_add_ps(euc1, euc2), _mm256_add_ps(euc3, euc4)));

for (size_t i = aligend_n; i < n; ++i)
{
const float num = p1[i] - p2[i];
distance += num * num * local_selected[i];
}

return distance;
}

Основной совет:

Не использовать uint16_t для вашего счетчика цикла, если вы действительно не хотите, чтобы компилятор урезал до 16 бит каждый раз. Используйте хотя бы unsigned, или вы иногда получаете лучше Asm от использования uintptr_t (или более условно, size_t). Нулевое расширение от 32-битной до ширины указателя происходит бесплатно на x86-64 только благодаря использованию 32-битных инструкций asm-размера с операндом, но иногда компиляторы все равно не справляются.

Используйте пять или более отдельных аккумуляторов вместо одного euclideanтаким образом, множество команд sub / FMA могут быть в полете без узких мест в задержке цепочки зависимостей, переносимых циклами, которая превращает FMA в один аккумулятор.

Время ожидания FMA в Intel Haswell составляет 5 циклов, а пропускная способность — один на 0,5 цикла. Смотрите также латентность против пропускной способности в Intel intrinsics, а также мой ответ на Почему Мулсс занимает всего 3 цикла в Haswell, в отличие от таблиц инструкций Агнера? для более продвинутой версии.

---

Избегайте передачи аргументов через глобальные переменные. Видимо твой n является константой времени компиляции (что хорошо), но selected_dimensions не так ли? Если это так, то вы используете только один набор масок во всей вашей программе, поэтому не забывайте о материалах, приведенных ниже о сжатии масок.

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

---

обновление: ваши массивы маленькие, всего ~ 100 элементов, так что развертывание только на 2 может быть хорошим, чтобы уменьшить накладные расходы при запуске / очистке. Внеочередное выполнение может скрыть задержку FMA на этом коротком числе итераций, особенно если конечный результат этого вызова функции не требуется для определения входных параметров для следующего вызова.

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

Как обсуждено в комментариях, очистка первой итерации цикла позволяет избежать первого FMA путем инициализации euc1 = stuff(p1[0], p2[0]); вместо _mm256_setzero_ps(),

Дополнение массивов к полному вектору (или даже к полному развернутому телу цикла из 2 векторов) с нулями позволяет полностью избежать скалярного цикла очистки и сделать всю функцию очень компактной.

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

Вы не показываете свой hsum256_ps_avx код, но это приличная часть общей задержки и, возможно, пропускной способности вашей функции. Убедитесь, что вы оптимизировали его для пропускной способности: например, избежать haddps / _mm_hadd_ps, Смотрите мой ответ на Самый быстрый способ сделать горизонтальную векторную сумму с плавающей точкой на x86.

---

Ваш конкретный случай:

Таким образом, я мог бы предварительно вычислить массив __m256 для передачи _mm256_blendv_ps вместо умножения на 0 или 1 в FMA.

Да, это было бы лучше, особенно если это позволяет вам свернуть что-то еще в FMAdd / FMSub. Но даже лучше, используйте логическое значение _mm256_and_ps со всеми нулями или со всеми единицами. Это оставляет значение без изменений (1 & x == x) или обнуляется (0 & x == 0и бинарное представление float 0.0 это все нули.)

Если ваши маски не пропали в кеше, храните их полностью распакованными, чтобы их можно было просто загрузить.

Если вы используете разные маски с одинаковыми p1 а также p2Вы могли бы предварительно вычислить p1-p2 в квадрате, а затем просто сделать в маске add_ps сокращение над этим. (Но обратите внимание, что FMA имеет лучшую пропускную способность, чем ADD на Intel pre-Skylake. У Haswell / Broadwell есть 2 модуля FMA, но запускаются ADDPS на выделенном модуле с меньшей задержкой (3c против 5c). Есть только один модуль добавления вектора-FP. Skylake просто запускает все на устройствах FMA с задержкой в ​​4 цикла.) В любом случае, это означает, что использование FMA в качестве 1.0 * x + y, Но вы, вероятно, в порядке, потому что вам все еще нужно загрузить маску и square(p1-p2) отдельно, так что это 2 загрузки на одно добавление FP, так что один за цикл не отстает от пропускной способности загрузки. Если только вы (или компилятор) не очистите несколько итераций впереди и сохраните плавающие данные для этих итераций в регистрах по нескольким различным local_selected маски.

---

ОбновитьЯ написал это, предполагая, что размер массива составляет 2-3 миллиона, а не ~ 100. Профиль для L1D-кэша не решает, стоит ли тратить больше инструкций ЦП для уменьшения объема кеш-памяти. Если вы всегда используете одну и ту же маску для всех 3 миллионов вызовов, вероятно, ее не стоит сжимать.

Вы можете сжать свои маски до 8 бит на элемент и загрузить их pmovsx (_mm256_cvtepi8_epi32) (Расширение знака «все единицы» дает более широкое значение «все единицы», потому что так дополняют 2 -1 работает). К сожалению, использование его в качестве нагрузки раздражает; компиляторы иногда не могут оптимизировать _mm256_cvtepi8_epi32(_mm_cvtsi64x_si128(foo)) в vpmovsxbd ymm0, [mem]и вместо этого на самом деле использовать отдельный vmovq инструкция.

const uint64_t *local_selected = something;  // packed to 1B per element

__m256 euc1 = _mm256_setzero_ps(), euc2 = _mm256_setzero_ps(),
euc3 =  _mm256_setzero_ps(), euc4 =  _mm256_setzero_ps();

for (i = 0 ; i < n ; i += 8*4) {  // 8 floats * an unroll of 4

__m256 mask = _mm256_castsi256_ps( _mm256_cvtepi8_epi32(_mm_cvtsi64x_si128(local_selected[i*1 + 0])) );
// __m256 mask = _mm256_load_ps(local_selected[i*8 + 0]); //  without packing

const __m256 r1 = _mm256_sub_ps(_mm256_load_ps(&p1[i*8 + 0]), _mm256_load_ps(&p2[i*8 + 0]));
r1 = _mm256_and_ps(r1, mask);             // zero r1 or leave it untouched.
euc1 = _mm256_fmadd_ps(r1, r1, euc1);    // euc1 += r1*r1
// ... same for r2 with local_selected[i + 1]
// and p1/p2[i*8 + 8]
// euc2 += (r2*r2) & mask2

// and again for euc3 (local_selected[i + 2], p1/p2[i*8 + 16]
// and again for euc3 (local_selected[i + 3], p1/p2[i*8 + 24]
}
euclidean = hsum (euc1+euc2+euc3+euc4);

Я предполагаю, что вы немного узкие места на пропускную способность без pmovsx, так как у вас есть три нагрузки для трех векторных операций ALU. (И с микро-слиянием, это всего 4 мопа слитых доменов на процессоре Intel, так что это не узкое место на переднем конце). И три ALU мопа могут работать на разных портах (vandps это 1 моп для порта 5 на Intel pre-Skylake. На SKL он может работать на любом порту).

Добавление в случайном порядке ( pmovsx) потенциально узкие места на порту 5 (на Haswell / Broadwell). Вы можете захотеть использовать vpand для маскировки, чтобы он мог работать на любом порту, если вы настраиваете на HSW / BDW, даже если они имеют дополнительную задержку обхода между математическими инструкциями целых чисел AND и FP. С достаточным количеством аккумуляторов вы не будете привязаны к задержке. (Skylake имеет дополнительную задержку обхода для VANDPS в зависимости от того, на каком порту он работает).

blendv медленнее AND: всегда как минимум 2 мопа.

---

Сжатие маски еще больше для больших массивов

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

Я думаю, что идеальный формат для ваших данных маски — это чередование 32 векторов масок, что делает его очень дешевым «распаковывать» на лету. Используйте сдвиг, чтобы ввести правильную маску в старший бит каждого 32-битного элемента, и используйте ее с vblendvps условно ноль элементов путем смешивания с нулем. (Или с арифметическим сдвигом вправо + логическое И)

__m256i masks = _mm256_load_si256(...);

// this actually needs a cast to __m256, omitted for readability
r0 = _mm256_blendv_ps(_mm256_setzero_ps(), r0, masks);
...

__m256i mask1 = _mm256_slli_epi32(masks, 1);
r1 = _mm256_blendv_ps(_mm256_setzero_ps(), r1, mask1);
...

__m256i mask2 = _mm256_slli_epi32(masks, 2);
r2 = _mm256_blendv_ps(_mm256_setzero_ps(), r2, mask2);
...

// fully unrolling is overkill; you can set up for a loop back to r0 with
masks = _mm256_slli_epi32(masks, 4);

Вы могли бы также сделать masks = _mm256_slli_epi32(masks, 1); на каждом шаге, что может быть лучше, потому что он использует 1 регистр меньше. Но он может быть более чувствительным к конфликтам ресурсов, вызывающим задержку в цепочке удаления маски, поскольку каждая маска зависит от предыдущей.

Intel Haswell работает как vblendvps моп только на port5, так что вы можете рассмотреть возможность использования _mm256_srai_epi32 + _mm256_and_ps, Но Skylake может запустить 2 мопа на любом из p015, так что смешивание здесь хорошо (хотя оно связывает векторный регистр, содержащий вектор со всеми нулями).

---

Генерация масок в этом чередованном формате с упакованным сравнением, затем _mm256_srli_epi32(cmp_result, 31) ИЛИ это в векторе, который вы строите. Затем сдвиньте его влево на 1. Повторите 32 раза.

---

Вы по-прежнему можете использовать этот формат, если в ваших массивах содержится менее 32 полных векторов данных. Младшие биты просто останутся неиспользованными. Или вы могли бы иметь маски для 2 или более selected_dimensions по вектору. например верхние 16 бит каждого элемента предназначены для одного selected_dimensionsи нижние 16 бит для другого. Вы могли бы сделать что-то вроде

__m256i masks =  _mm256_load_si256(dimensions[selector/2]);
masks = _mm256_sll_epi32(masks, 16 * (selector % 2));

// or maybe
if (selector % 2) {
masks = _mm256_slli_epi32(masks, 16);
}

---

AVX512:

AVX512 может напрямую использовать растровую маску, поэтому она несколько более эффективна. Просто используйте const __mmask16 *local_selected = whatever; объявить массив 16-битных масок (для использования с 512b векторами из 16 чисел с плавающей запятой) и использовать r0 = _mm512_maskz_sub_ps(p1,p2, local_selected[i]); чтобы обнулить вычитание.

Если вы на самом деле являетесь узким местом в пропускной способности uop загрузочного порта (2 загрузки в такт), вы можете попробовать загрузить 64 бита данных маски одновременно и использовать сдвиг маски, чтобы получить другой низкий 16 из них. Это, вероятно, не будет проблемой, если ваши данные не будут загружены в кэш L1D.

Во-первых, очень просто сгенерировать данные маски с помощью сравнения в маску без чередования.

---

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