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

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

本文关键字：矢量化 更好 性能 边缘检测 何进行 代码      更新时间：2024-05-23

我一直在学习这门课程，在某个时候给出了下面的代码，讲师声称矢量化是通过在内部和外部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

这里有一些未经测试的概念验证实现，它每个数据包使用4个加法、1个fmsub和3个加载(而不是直接实现的9个加载、7个加法和1个fmsub(。我省略了箝位(至少对于float图像来说，这看起来很不寻常，对于uint8来说，它什么都没做，除非你把P val = ...改成auto val = ...，正如Peter在评论中注意到的那样(——但你可以很容易地自己添加它。

该实现的思想是将左右像素(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。。。像素以获得更好的缓存位置(只要您有足够的寄存器(。当然，这会使实现这两个功能稍微复杂一些
  • 如果展开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实现与编译器从通用实现生成的内容。

相关文章：