Как выполнить векторизацию кода вручную с более высокой производительностью, чем автоматическая c векторизация для обнаружения краев - PullRequest
Новая
       31

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

2 голосов
/ 07 мая 2020
• 1000 трудно. Как я могу самостоятельно векторизовать код, используемый в курсе, и есть ли способ добиться большей производительности, чем если бы я просто добавил #pragma omp simd и двинулся дальше? 
template<typename P>
void ApplyStencil(ImageClass<P> & img_in, ImageClass<P> & img_out) {

  const int width  = img_in.width;
  const int height = img_in.height;

  P * in  = img_in.pixel;
  P * out = img_out.pixel;

  for (int i = 1; i < height-1; i++)
    for (int j = 1; j < width-1; j++) {
      P val = -in[(i-1)*width + j-1] -   in[(i-1)*width + j] - in[(i-1)*width + j+1] 
    -in[(i  )*width + j-1] + 8*in[(i  )*width + j] - in[(i  )*width + j+1] 
    -in[(i+1)*width + j-1] -   in[(i+1)*width + j] - in[(i+1)*width + j+1];

      val = (val < 0   ? 0   : val);
      val = (val > 255 ? 255 : val);

      out[i*width + j] = val;
    }

}

template void ApplyStencil<float>(ImageClass<float> & img_in, ImageClass<float> & img_out);

Я компилирую, используя g cc с флагами -march=native -fopenmp для поддержки AVX512 на процессоре Skylake. 

❯ gcc -march=native -Q --help=target|grep march
  -march=                           skylake

❯ gcc -march=knl -dM -E - < /dev/null | egrep "SSE|AVX" | sort
#define __AVX__ 1
#define __AVX2__ 1
#define __AVX512CD__ 1
#define __AVX512ER__ 1
#define __AVX512F__ 1
#define __AVX512PF__ 1
#define __SSE__ 1
#define __SSE2__ 1
#define __SSE2_MATH__ 1
#define __SSE3__ 1
#define __SSE4_1__ 1
#define __SSE4_2__ 1
#define __SSE_MATH__ 1
#define __SSSE3__ 1

1 Ответ

2 голосов
/ 09 мая 2020

Вот некоторая непроверенная реализация, подтверждающая концепцию, которая использует 4 добавления, 1 fmsub и 3 загрузки на пакет (вместо 9 загрузок, 7 добавлений, 1 fmsub для прямой реализации). Я не учел зажим (который для изображений float выглядит по крайней мере необычно, а для uint8 он ничего не делает, если вы не измените P val = ... на auto val = ..., как заметил Питер в комментариях), но вы можете легко добавьте это самостоятельно.

Идея этой реализации состоит в том, чтобы суммировать пиксели слева и справа (x0_2), а также все 3 (x012) и складывать их из 3 последовательных строк (a012 + b0_2 + c012 ), затем вычтите это из среднего пикселя, умноженного на 8. В конце каждого l oop опустите содержимое a012 и переместите bX в aX и cX в bX для следующей итерации.

Функция applyStencil просто применяет первую функцию для каждого столбца из 16 пикселей (начиная с col = 1 и в конце просто выполняет возможное перекрытие вычислений для последних 16 столбцов). Если ваше входное изображение имеет менее 18 столбцов, вам нужно обработать это по-другому (возможно, с помощью маскированных загрузок / сохранений). 

#include <immintrin.h>

void applyStencilColumn(float const *in, float *out, size_t width, size_t height)
{
  if(height < 3) return; // sanity check
  float const* last_in = in + height*width;
  __m512 a012, b012, b0_2, b1;
  __m512 const eight = _mm512_set1_ps(8.0);
  {
    // initialize first rows:
    __m512 a0 = _mm512_loadu_ps(in-1);
    __m512 a1 = _mm512_loadu_ps(in+0);
    __m512 a2 = _mm512_loadu_ps(in+1);
    a012 = _mm512_add_ps(_mm512_add_ps(a0,a2),a1);
    in += width;
    __m512 b0 = _mm512_loadu_ps(in-1);
    b1 = _mm512_loadu_ps(in+0);
    __m512 b2 = _mm512_loadu_ps(in+1);
    b0_2 = _mm512_add_ps(b0,b2);
    b012 = _mm512_add_ps(b0_2,b1);
    in += width;
  }
  // skip first row for output:
  out += width;

  for(; in<last_in; in+=width, out+=width)
  {
    // precalculate sums for next row:
    __m512 c0 = _mm512_loadu_ps(in-1);
    __m512 c1 = _mm512_loadu_ps(in+0);
    __m512 c2 = _mm512_loadu_ps(in+1);
    __m512 c0_2 = _mm512_add_ps(c0,c2);
    __m512 c012 = _mm512_add_ps(c0_2, c1);

    __m512 outer = _mm512_add_ps(_mm512_add_ps(a012,b0_2), c012);
    __m512 result = _mm512_fmsub_ps(eight, b1, outer);

    _mm512_storeu_ps(out, result);
    // shift/rename registers (with some unrolling this can be avoided entirely)
    a012 = b012;
    b0_2 = c0_2; b012 = c012; b1 = c1;
  }
}

void applyStencil(float const *in, float *out, size_t width, size_t height)
{
  if(width < 18) return; // assert("special case of narrow image not implemented");

  for(size_t col = 1; col < width - 18; col += 16)
  {
    applyStencilColumn(in + col, out + col, width, height);
  }
  applyStencilColumn(in + width - 18, out + width - 18, width, height);
}

Возможные улучшения (оставлены в качестве упражнения):

  • 1025 * может воздействовать на столбцы размером 32, 48, 64, ... пикселей для лучшей локализации кеша (если у вас достаточно регистров). Это, конечно, немного усложняет реализацию обеих функций.
  • Если вы развернете 3 (или 6, 9, ...) итерации for(; in<last_in; in+=width) l oop, в действительности не потребуется перемещать регистры (плюс общее преимущество разворачивания).
  • Если ваша ширина кратна 16, вы можете убедиться, что по крайней мере хранилища в основном выровнены (за исключением первого и последнего столбцов).
  • Вы можете выполнять итерацию по небольшому количеству строк одновременно (добавив еще один внешний l oop к основной функции и вызвав applyStencilColumn с фиксированной высотой. Убедитесь, что у вас есть правильная степень перекрытия между наборы строк. (Идеальное количество строк, вероятно, зависит от размера вашего изображения).
  • Вы также всегда можете добавить 3 последовательных пикселя, но вместо этого умножьте центральный пиксель на 9 (9*b1-outer). Затем ( с некоторыми бухгалтерскими усилиями) вы можете добавить row0+(row1+row2) и (row1+row2)+row3, чтобы получить промежуточные результаты row1 и row2 (имея 3 прибавления вместо 4). То же самое по горизонтали выглядит однако более сложный.

Конечно, вы всегда должны тестировать и тестировать любую пользовательскую реализацию SIMD по сравнению с тем, что ваш компилятор генерирует из общей реализации c.

...