如何以比边缘检测自动矢量化更好的性能进行手动代码矢量化

Question

我一直在学习 this coursera 课程，在某些时候给出了下面的代码，讲师声称通过在内部和外部 for 之间包含 #pragma omp simd 来完成矢量化循环，因为引导向量化很难。我如何自己对课程中使用的代码进行向量化处理，有没有一种方法可以比我简单地添加 #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);

我正在使用 gcc 和 -march=native -fopenmp 标志进行编译，以便在 skylake 处理器上支持 AVX512。

❯ 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

Answer 1

这是一些未经测试的概念验证实现，每个数据包使用 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 列，您需要以不同的方式处理（可能通过屏蔽 loads/stores）。

#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、... 像素的列以获得更好的缓存局部性（只要您有足够的寄存器）。当然，这会使这两个功能的实现稍微复杂一些。
如果展开 for(; in<last_in; in+=width) 循环的 3（或 6、9、...）次迭代，则无需实际移动寄存器（加上展开的一般好处）。
如果您的宽度是 16 的倍数，您可以确保至少商店大部分对齐（第一列和最后一列除外）。
您可以同时迭代少量行（通过向主函数添加另一个外部循环并使用固定高度调用 applyStencilColumn。确保有适当的重叠量行集之间。（理想的行数可能取决于图像的大小）。
您也可以始终添加 3 个连续像素，但将中心像素乘以 9 (9*b1-outer)。然后（通过一些簿记工作）您可以添加 row0+(row1+row2) 和 (row1+row2)+row3 以获得 row1 和 row2 中间结果（有 3 个而不是 4 个加法）。不过，水平地做同样的事情看起来更复杂。

当然，您应该始终对任何自定义 SIMD 实现与您的编译器从通用实现生成的实现进行测试和基准测试。

如何以比边缘检测自动矢量化更好的性能进行手动代码矢量化

How to do manual code vectorization with better performance that automatic vectorization for edge detection

c++

optimization

avx512