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

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

Идея этой реализации состоит в том, чтобы суммировать пиксели слева и справа (x0_2), а также все 3 (x012) и сложить их из 3 последовательных строк (a012 + b0_2 + c012), а затем вычесть это из среднего пикселя, умноженного на 8. В конце каждый цикл удаляет содержимое 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);
}

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

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

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