SIMD для работы с плавающим порогом - PullRequest
Новая
       50

SIMD для работы с плавающим порогом

3 голосов
/ 19 июня 2019

Я хотел бы ускорить вычисление векторов, и я считаю, что SIMD-инструкции для сравнения и манипулирования числами могут помочь, вот операция: 

void func(const double* left, const double* right, double* res, const size_t size, const double th, const double drop) {
        for (size_t i = 0; i < size; ++i) {
            res[i] = right[i] >= th ? left[i] : (left[i] - drop) ;
        }
    }

В основном, это сбрасывает значение leftна drop в случае, если right значение больше, чем threshold.

Размер составляет около 128-256 (не такой большой), но вычисления вызываются интенсивно.

Я пыталсяначать с развертывания цикла, но не выиграл много производительности, но, возможно, необходимы некоторые инструкции по компиляции.

Не могли бы вы предложить некоторые улучшения в коде для более быстрого вычисления?

Ответы [ 3 ]

5 голосов
/ 20 июня 2019

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

К сожалению, gcc только автоматически векторизуется с -ffast-math. Я не уверен, какая часть строгого FP останавливает это. (Обновление: требуется только -fno-trapping-math: это, вероятно, безопасно в большинстве случаев, если вы не снимаете маску с каких-либо исключений FP или не смотрите на флаги исключений MXCSR sticky FP.)

С этой опцией GCC также будет использовать (v)pblendvpd с -march=nehalem или -march=znver1. Посмотри на Годболт

Кроме того, ваша функция C не работает. th и drop являются скалярными двойными, но вы объявляете их как const double *

AVX512F позволит вам сравнить !(right[i] >= thresh) и использовать полученную маску для вычитания с маскированием слиянием.

Элементы, где предикат был истинным, получат left[i] - drop, другие элементы сохранят свое значение left[i], потому что вы объединяете информацию в вектор left значений.

К сожалению, GCC с -march=skylake-avx512 использует обычный vsubpd, а затем отдельный vmovapd zmm2{k1}, zmm5 для смешивания, что, очевидно, является пропущенной оптимизацией. Назначение наложения уже является одним из входов SUB.

Использование AVX512VL для 256-битных векторов (в случае, если остальная часть вашей программы не может эффективно использовать 512-битный, поэтому вы не страдаете от снижения турбо тактовой частоты): 

__m256d left = ...;
__m256d right = ...;
__mmask8 cmp = _mm256_cmp_pd_mask(right, set1(th), _CMP_NGE_UQ);
__m256d res = _mm256_mask_sub_pd (left, cmp, left, set1(drop));

Итак (помимо загрузки и хранения) это 2 инструкции с AVX512F / VL.

Если вам не нужно определенное поведение NaN в вашей версии, GCC также может автоматически векторизовать

И это более эффективно для всех компиляторов, потому что вам просто нужно AND, а не переменное смешение. Так что это значительно лучше только с SSE2, а также с большинством процессоров, даже если они поддерживают SSE4. 1 blendvpd, потому что эта инструкция не так эффективна.

Вы можете вычесть 0.0 или drop из left[i] на основе результата сравнения.

Создание 0.0 или константы на основе результата сравнения чрезвычайно эффективно: просто инструкция andps. (Битовый шаблон для 0.0 - это все нули, и SIMD сравнивает производящие векторы со всеми 1 или 0 битами. Таким образом, И сохраняет старое значение или обнуляет его.)

Мы также можем добавить -drop вместо вычитания drop. Это требует дополнительного отрицания при вводе, но с AVX позволяет операнду источника памяти для vaddpd. GCC предпочитает использовать режим индексированной адресации, чтобы на самом деле не уменьшить количество внешних операций на процессорах Intel; это будет "не ламинировать". Но даже с -ffast-math gcc не выполняет эту оптимизацию самостоятельно, чтобы разрешить сворачивание нагрузки. (Хотя не стоит делать отдельные приращения указателя, если мы не развернем цикл.) 

void func3(const double *__restrict left, const double *__restrict right, double *__restrict res,
  const size_t size, const double th, const double drop)
{
    for (size_t i = 0; i < size; ++i) {
        double add = right[i] >= th ? 0.0 : -drop;
        res[i] = left[i] + add;
    }
}
Внутренний цикл

GCC 9.1 (без каких-либо опций -march и без -ffast-math) из ссылки Godbolt выше: 

# func3 main loop
# gcc -O3 -march=skylake       (without fast-math)
.L33:
    vcmplepd        ymm2, ymm4, YMMWORD PTR [rsi+rax]
    vandnpd ymm2, ymm2, ymm3
    vaddpd  ymm2, ymm2, YMMWORD PTR [rdi+rax]
    vmovupd YMMWORD PTR [rdx+rax], ymm2
    add     rax, 32
    cmp     r8, rax
    jne     .L33

Или простая версия SSE2 имеет внутренний цикл, который в основном такой же, как с left - zero_or_drop вместо left + zero_or_minus_drop, поэтому, если вы не можете пообещать компилятору 16-байтовое выравнивание или вы делаете версию AVX, отрицая drop это просто дополнительные накладные расходы.

Отрицание drop берет константу из памяти (для XOR знакового бита), и это единственная статическая константа, которая нужна этой функции , так что компромисс стоит рассмотреть в вашем случае, когда цикл не ' Я бегаю огромное количество раз. (Если th или drop также не являются константами времени компиляции после вставки и все равно загружаются. Или особенно если -drop можно вычислить во время компиляции. Или если вы можете заставить вашу программу работать с отрицательным drop.)

Другое различие между сложением и вычитанием состоит в том, что вычитание не уничтожает знак нуля. -0.0 - 0.0 = -0.0, +0.0 - 0.0 = +0.0. В случае, если это имеет значение. 

# gcc9.1 -O3
.L26:
    movupd  xmm5, XMMWORD PTR [rsi+rax]
    movapd  xmm2, xmm4                    # duplicate  th
    movupd  xmm6, XMMWORD PTR [rdi+rax]
    cmplepd xmm2, xmm5                    # destroy the copy of th
    andnpd  xmm2, xmm3                    # _mm_andnot_pd
    addpd   xmm2, xmm6                    # _mm_add_pd
    movups  XMMWORD PTR [rdx+rax], xmm2
    add     rax, 16
    cmp     r8, rax
    jne     .L26

GCC использует невыровненные нагрузки, поэтому (без AVX) он не может сложить операнд источника памяти в cmppd или subpd

4 голосов
/ 20 июня 2019

Вот, пожалуйста (непроверенный), я пытался объяснить в комментариях, что они делают. 

void func_sse41( const double* left, const double* right, double* res,
    const size_t size, double th, double drop )
{
    // Verify the size is even.
    // If it's not, you'll need extra code at the end to process last value the old way.
    assert( 0 == ( size % 2 ) );

    // Load scalar values into 2 registers.
    const __m128d threshold = _mm_set1_pd( th );
    const __m128d dropVec = _mm_set1_pd( drop );

    for( size_t i = 0; i < size; i += 2 )
    {
        // Load 4 double values into registers, 2 from right, 2 from left
        const __m128d r = _mm_loadu_pd( right + i );
        const __m128d l = _mm_loadu_pd( left + i );
        // Compare ( r >= threshold ) for 2 values at once
        const __m128d comp = _mm_cmpge_pd( r, threshold );
        // Compute ( left[ i ] - drop ), for 2 values at once
        const __m128d dropped = _mm_sub_pd( l, dropVec );
        // Select either left or ( left - drop ) based on the comparison.
        // This is the only instruction here that requires SSE 4.1.
        const __m128d result = _mm_blendv_pd( l, dropped, comp );
        // Store the 2 result values
        _mm_storeu_pd( res, result );
    }
}

Код завершится с ошибкой «недопустимая инструкция», если у ЦПУ нет SSE 4.1. Для достижения наилучшего результата, обнаружение с идентификатором процессора изящно сбои. Я думаю, что сейчас, в 2019 году, вполне разумно предположить, что эта поддержка поддерживается, Intel сделала в 2008 году, AMD - в 2011 году, согласно исследованию steam, «96,3%». Если вы хотите поддерживать старые процессоры, можно эмулировать _mm_blendv_pd с помощью 3 других инструкций: _mm_and_pd, _mm_andnot_pd, _mm_or_pd.

Если вы можете гарантировать, что данные выровнены, замена нагрузок на _mm_load_pd будет немного быстрее, _mm_cmpge_pd компилируется в CMPPD https://www.felixcloutier.com/x86/cmppd, который может принимать один из аргументов непосредственно из ОЗУ.

Потенциально, вы можете получить дальнейшее 2-х кратное улучшение, написав версию AVX. Но я надеюсь, что даже версия SSE быстрее, чем ваш код, она обрабатывает 2 значения за итерацию и не имеет условий внутри цикла. Если вам не повезло, AVX будет работать медленнее, многим процессорам потребуется некоторое время для включения их устройств AVX, что займет много тысяч циклов. До включения код AVX работает очень медленно.

2 голосов
/ 20 июня 2019

Вы можете использовать векторные расширения GCC и Clang для реализации троичной функции выбора (см. https://stackoverflow.com/a/48538557/2542702).

#include <stddef.h>
#include <inttypes.h>

#if defined(__clang__)
typedef  double double4 __attribute__ ((ext_vector_type(4)));
typedef int64_t   long4 __attribute__ ((ext_vector_type(4)));
#else
typedef  double double4 __attribute__ ((vector_size (sizeof(double)*4)));
typedef int64_t   long4 __attribute__ ((vector_size (sizeof(int64_t)*4)));
#endif

double4 select(long4 s, double4 a, double4 b) {
  double4 c;
  #if defined(__GNUC__) && !defined(__INTEL_COMPILER) && !defined(__clang__)
  c = s ? a : b;
  #else
  for(int i=0; i<4; i++) c[i] = s[i] ? a[i] : b[i];
  #endif
  return c;
}

void func(double* left, double* right, double* res, size_t size, double th, double drop) {
  size_t i;
  for (i = 0; i<(size&-4); i+=4) {
    double4 leftv = *(double4*)&left[i];
    double4 rightv = *(double4*)&right[i];
    *(double4*)&res[i] = select(rightv >= th, leftv, leftv - drop);
  }
  for(;i<size; i++) res[i] = right[i] >= th ? left[i] : (left[i] - drop);
}

https://godbolt.org/z/h4OKMl

...