bgp*_*000 2 c++ x86 sse image-processing simd
我有一个功能可以将 8 位图像缩小两倍。我之前已经用 SSE 优化了 rgb32 案例。现在我想对 gray8 案例做同样的事情。
在核心，有一个函数取两行像素数据，其工作方式如下：
/**
* Calculates the average of two rows of gray8 pixels by averaging four pixels.
*/
void average2Rows(const uint8_t* row1, const uint8_t* row2, uint8_t* dst, int size)
{
for (int i = 0; i < size - 1; i += 2)
*(dst++) = ((row1[i]+row1[i+1]+row2[i]+row2[i+1])/4)&0xFF;
}
现在，我想出了一个 SSE 变体，它大约快三倍，但它确实涉及很多改组，我认为可能会做得更好。有人看到这里可以优化什么吗？
/* row1: 16 8-bit values A-P
* row2: 16 8-bit values a-p
* returns 16 8-bit values (A+B+a+b)/4, (C+D+c+d)/4, ..., (O+P+o+p)/4
*/
__m128i avg16Bytes(const __m128i& row1, const __m128i& row2)
{
static const __m128i zero = _mm_setzero_si128();
__m128i ABCDEFGHIJKLMNOP = _mm_avg_epu8(row1_u8, row2);
__m128i ABCDEFGH = _mm_unpacklo_epi8(ABCDEFGHIJKLMNOP, zero);
__m128i IJKLMNOP = _mm_unpackhi_epi8(ABCDEFGHIJKLMNOP, zero);
__m128i AIBJCKDL = _mm_unpacklo_epi16( ABCDEFGH, IJKLMNOP );
__m128i EMFNGOHP = _mm_unpackhi_epi16( ABCDEFGH, IJKLMNOP );
__m128i AEIMBFJN = _mm_unpacklo_epi16( AIBJCKDL, EMFNGOHP );
__m128i CGKODHLP = _mm_unpackhi_epi16( AIBJCKDL, EMFNGOHP );
__m128i ACEGIKMO = _mm_unpacklo_epi16( AEIMBFJN, CGKODHLP );
__m128i BDFHJLNP = _mm_unpackhi_epi16( AEIMBFJN, CGKODHLP );
return _mm_avg_epu8(ACEGIKMO, BDFHJLNP);
}
/*
* Calculates the average of two rows of gray8 pixels by averaging four pixels.
*/
void average2Rows(const uint8_t* src1, const uint8_t* src2, uint8_t* dst, int size)
{
for(int i = 0;i<size-31; i+=32)
{
__m128i tl = _mm_loadu_si128((__m128i const*)(src1+i));
__m128i tr = _mm_loadu_si128((__m128i const*)(src1+i+16));
__m128i bl = _mm_loadu_si128((__m128i const*)(src2+i));
__m128i br = _mm_loadu_si128((__m128i const*)(src2+i+16)))
__m128i l_avg = avg16Bytes(tl, bl);
__m128i r_avg = avg16Bytes(tr, br);
_mm_storeu_si128((__m128i *)(dst+(i/2)), _mm_packus_epi16(l_avg, r_avg));
}
}
笔记：
编辑：现在有一个github 存储库实现了这个问题的答案。最快的解决方案是由用户Peter Cordes提供的。有关详细信息，请参阅下面的他的文章：
__m128i avg16Bytes(const __m128i& row1, const __m128i& row2)
{
// Average the first 16 values of src1 and src2:
__m128i avg = _mm_avg_epu8(row1, row2);
// Unpack and horizontal add:
avg = _mm_maddubs_epi16(avg, _mm_set1_epi8(1));
// Divide by 2:
return _mm_srli_epi16(avg, 1);
}
它通过计算(a+b)/2 + (c+d)/2而不是 来作为我的原始实现(a+b+c+d)/4，因此它具有相同的一对一舍入误差。
感谢用户Paul R实施了一个比我快两倍但准确的解决方案：
__m128i avg16Bytes(const __m128i& row1, const __m128i& row2)
{
// Unpack and horizontal add:
__m128i row1 = _mm_maddubs_epi16(row1_u8, _mm_set1_epi8(1));
__m128i row2 = _mm_maddubs_epi16(row2_u8, _mm_set1_epi8(1));
// vertical add:
__m128i avg = _mm_add_epi16(row1_avg, row2_avg);
// divide by 4:
return _mm_srli_epi16(avg, 2);
}
如果您愿意接受pavgb两次使用的双舍入，您可以通过首先使用 进行垂直平均pavgb，将需要解压缩为 16 位元素的数据量减少一半，从而比 Paul R 的答案更快。（并允许一半的负载折叠到内存操作数中pavgb，从而减少某些 CPU 的前端瓶颈。）
对于水平平均，最好的办法是可能仍然pmaddubsw与set1(1)和1移位，然后收拾。
// SSSE3 version
// I used `__restrict__` to give the compiler more flexibility in unrolling
void average2Rows_doubleround(const uint8_t* __restrict__ src1, const uint8_t*__restrict__ src2,
uint8_t*__restrict__ dst, size_t size)
{
const __m128i vk1 = _mm_set1_epi8(1);
size_t dstsize = size/2;
for (size_t i = 0; i < dstsize - 15; i += 16)
{
__m128i v0 = _mm_load_si128((const __m128i *)&src1[i*2]);
__m128i v1 = _mm_load_si128((const __m128i *)&src1[i*2 + 16]);
__m128i v2 = _mm_load_si128((const __m128i *)&src2[i*2]);
__m128i v3 = _mm_load_si128((const __m128i *)&src2[i*2 + 16]);
__m128i left = _mm_avg_epu8(v0, v2);
__m128i right = _mm_avg_epu8(v1, v3);
__m128i w0 = _mm_maddubs_epi16(left, vk1); // unpack and horizontal add
__m128i w1 = _mm_maddubs_epi16(right, vk1);
w0 = _mm_srli_epi16(w0, 1); // divide by 2
w1 = _mm_srli_epi16(w1, 1);
w0 = _mm_packus_epi16(w0, w1); // pack
_mm_storeu_si128((__m128i *)&dst[i], w0);
}
}
另一种选择是_mm_srli_epi16(v, 8)将奇数元素与每个水平对的偶数元素对齐。但是由于没有带有截断的水平包装，因此您必须_mm_and_si128(v, _mm_set1_epi16(0x00FF))在包装之前将其分成两半。事实证明，它比使用 SSSE3 慢pmaddubsw，尤其是在没有 AVX 的情况下，它需要额外的 MOVDQA 指令来复制寄存器。
void average2Rows_doubleround_SSE2(const uint8_t* __restrict__ src1, const uint8_t* __restrict__ src2, uint8_t* __restrict__ dst, size_t size)
{
size /= 2;
for (size_t i = 0; i < size - 15; i += 16)
{
__m128i v0 = _mm_load_si128((__m128i *)&src1[i*2]);
__m128i v1 = _mm_load_si128((__m128i *)&src1[i*2 + 16]);
__m128i v2 = _mm_load_si128((__m128i *)&src2[i*2]);
__m128i v3 = _mm_load_si128((__m128i *)&src2[i*2 + 16]);
__m128i left = _mm_avg_epu8(v0, v2);
__m128i right = _mm_avg_epu8(v1, v3);
__m128i l_odd = _mm_srli_epi16(left, 8); // line up horizontal pairs
__m128i r_odd = _mm_srli_epi16(right, 8);
__m128i l_avg = _mm_avg_epu8(left, l_odd); // leaves garbage in the high halves
__m128i r_avg = _mm_avg_epu8(right, r_odd);
l_avg = _mm_and_si128(l_avg, _mm_set1_epi16(0x00FF));
r_avg = _mm_and_si128(r_avg, _mm_set1_epi16(0x00FF));
__m128i avg = _mm_packus_epi16(l_avg, r_avg); // pack
_mm_storeu_si128((__m128i *)&dst[i], avg);
}
}
对于 AVX512BW，有_mm_cvtepi16_epi8，但 IACA 说它在 Skylake-AVX512 上是 2 uop，它只需要 1 个输入并产生半宽输出。根据 IACA，内存目标形式是总共 4 个未融合域 uops（与 reg,reg +单独存储相同）。我不得不使用_mm_mask_cvtepi16_storeu_epi8(&dst\[i+0\], -1, l_avg);它来获取它，因为 gcc 和 clang 无法将单独_mm_store的vpmovwb. （没有非屏蔽存储内在，因为编译器应该为你做这件事，就像他们_mm_load为典型的 ALU 指令折叠成内存操作数一样）。
它可能只在缩小到 1/4 或 1/8th ( cvtepi64_epi8)时有用，而不仅仅是缩小一半。或者可能有助于避免需要第二次 shuffle 来处理_mm512_packus_epi16. 使用 AVX2，在_mm256_packus_epi16on 之后[D C] [B A]，您有[D B | C A]，您可以使用 AVX2 修复它_mm256_permute4x64_epi64 (__m256i a, const int imm8)以在 64 位块中随机播放。但是对于 AVX512，你需要一个矢量 shuffle-control 来控制vpermq. packus+ 不过，修复洗牌可能仍然是更好的选择。
一旦你这样做了，循环中就没有多少向量指令了，而且让编译器使 asm 变得更紧会有很多好处。不幸的是，编译器很难做好你的循环。（这也有助于 Paul R 的解决方案，因为他从问题中复制了对编译器不友好的循环结构。）
以 gcc/clang 可以更好地优化的方式使用循环计数器，并使用避免每次通过循环重新进行符号扩展的类型。
对于您当前的循环，gcc/clang 实际上对 进行算术右移i/2，而不是增加 16（而不是 32）并使用缩放索引寻址模式进行加载。似乎他们没有意识到这i总是均匀的。
（Matt Godbolt 的编译器浏览器上的完整代码 + asm）：
.LBB1_2: ## clang's inner loop for int i, dst[i/2] version
movdqu xmm1, xmmword ptr [rdi + rcx]
movdqu xmm2, xmmword ptr [rdi + rcx + 16]
movdqu xmm3, xmmword ptr [rsi + rcx]
movdqu xmm4, xmmword ptr [rsi + rcx + 16]
pavgb xmm3, xmm1
pavgb xmm4, xmm2
pmaddubsw xmm3, xmm0
pmaddubsw xmm4, xmm0
psrlw xmm3, 1
psrlw xmm4, 1
packuswb xmm3, xmm4
mov eax, ecx # This whole block is wasted instructions!!!
shr eax, 31
add eax, ecx
sar eax # eax = ecx/2, with correct rounding even for negative `i`
cdqe # sign-extend EAX into RAX
movdqu xmmword ptr [rdx + rax], xmm3
add rcx, 32 # i += 32
cmp rcx, r8
jl .LBB1_2 # }while(i < size-31)
gcc7.1 并不是那么糟糕，（只是mov/ sar/ movsx），但是 gcc5.x 和 6.x 为 src1 和 src2 以及商店的计数器/索引做单独的指针增量。（完全是脑残行为，特别是因为他们仍然使用-march=sandybridge. 索引movdqu存储和非索引movdqu加载为您提供最大的循环开销。）
无论如何，在循环内使用dstsize和乘法i而不是除法会产生更好的结果。不同版本的 gcc 和 clang 可靠地将其编译为单个循环计数器，它们与负载的缩放索引寻址模式一起使用。你得到如下代码：
movdqa xmm1, xmmword ptr [rdi + 2*rax]
movdqa xmm2, xmmword ptr [rdi + 2*rax + 16]
pavgb xmm1, xmmword ptr [rsi + 2*rax]
pavgb xmm2, xmmword ptr [rsi + 2*rax + 16] # saving instructions with aligned loads, see below
...
movdqu xmmword ptr [rdx + rax], xmm1
add rax, 16
cmp rax, rcx
jb .LBB0_2
我曾经size_t i匹配size_t大小，以确保 gcc 不会浪费任何指令将符号扩展或零扩展到指针的宽度。（零扩展通常发生于免费的，虽然如此，unsigned size而且unsigned i可能已经确定，并保存一对夫妇REX前缀。）
您仍然可以摆脱cmp但将索引向上计算为 0，这将比我所做的更快地加快循环速度。我不确定让编译器不愚蠢并忽略cmp指令（如果您确实数到零）有多么容易。不过，从对象末尾开始索引没有问题。src1+=size;. 但是，如果您想使用未对齐的清理循环，它确实会使事情复杂化。
在我的 Skylake i7-6700k（最大 turbo 4.4GHz，但查看时钟周期计数而不是时间）。使用 g++7.1，对于 1024 字节的 100M 代表，这与 ~3.3 秒相差 ~2.7 秒。
Performance counter stats for './grayscale-dowscale-by-2.inline.gcc-skylake-noavx' (2 runs):
2731.607950 task-clock (msec) # 1.000 CPUs utilized ( +- 0.40% )
2 context-switches # 0.001 K/sec ( +- 20.00% )
0 cpu-migrations # 0.000 K/sec
88 page-faults:u # 0.032 K/sec ( +- 0.57% )
11,917,723,707 cycles # 4.363 GHz ( +- 0.07% )
42,006,654,015 instructions # 3.52 insn per cycle ( +- 0.00% )
41,908,837,143 uops_issued_any # 15342.186 M/sec ( +- 0.00% )
49,409,631,052 uops_executed_thread # 18088.112 M/sec ( +- 0.00% )
3,301,193,901 branches # 1208.517 M/sec ( +- 0.00% )
100,013,629 branch-misses # 3.03% of all branches ( +- 0.01% )
2.731715466 seconds time elapsed ( +- 0.40% )
与相同的矢量化，但具有int i并dst[i/2]创建更高的循环开销（更多标量指令）：
Performance counter stats for './grayscale-dowscale-by-2.loopoverhead-aligned-inline.gcc-skylake-noavx' (2 runs):
3314.335833 task-clock (msec) # 1.000 CPUs utilized ( +- 0.02% )
4 context-switches # 0.001 K/sec ( +- 14.29% )
0 cpu-migrations # 0.000 K/sec
88 page-faults:u # 0.026 K/sec ( +- 0.57% )
14,531,925,552 cycles # 4.385 GHz ( +- 0.06% )
51,607,478,414 instructions # 3.55 insn per cycle ( +- 0.00% )
51,109,303,460 uops_issued_any # 15420.677 M/sec ( +- 0.00% )
55,810,234,508 uops_executed_thread # 16839.040 M/sec ( +- 0.00% )
3,301,344,602 branches # 996.080 M/sec ( +- 0.00% )
100,025,451 branch-misses # 3.03% of all branches ( +- 0.00% )
3.314418952 seconds time elapsed ( +- 0.02% )
与 Paul R 的版本（针对较低的循环开销进行了优化）：准确但速度较慢
Performance counter stats for './grayscale-dowscale-by-2.paulr-inline.gcc-skylake-noavx' (2 runs):
3751.990587 task-clock (msec) # 1.000 CPUs utilized ( +- 0.03% )
3 context-switches # 0.001 K/sec
0 cpu-migrations # 0.000 K/sec
88 page-faults:u # 0.024 K/sec ( +- 0.56% )
16,323,525,446 cycles # 4.351 GHz ( +- 0.04% )
58,008,101,634 instructions # 3.55 insn per cycle ( +- 0.00% )
57,610,721,806 uops_issued_any # 15354.709 M/sec ( +- 0.00% )
55,505,321,456 uops_executed_thread # 14793.566 M/sec ( +- 0.00% )
3,301,456,435 branches # 879.921 M/sec ( +- 0.00% )
100,001,954 branch-misses # 3.03% of all branches ( +- 0.02% )
3.752086635 seconds time elapsed ( +- 0.03% )
与带有额外循环开销的Paul R 的原始版本对比：
Performance counter stats for './grayscale-dowscale-by-2.loopoverhead-paulr-inline.gcc-skylake-noavx' (2 runs):
4154.300887 task-clock (msec) # 1.000 CPUs utilized ( +- 0.01% )
3 context-switches # 0.001 K/sec
0 cpu-migrations # 0.000 K/sec
90 page-faults:u # 0.022 K/sec ( +- 1.68% )
18,174,791,383 cycles # 4.375 GHz ( +- 0.03% )
67,608,724,157 instructions # 3.72 insn per cycle ( +- 0.00% )
66,937,292,129 uops_issued_any # 16112.769 M/sec ( +- 0.00% )
61,875,610,759 uops_executed_thread # 14894.350 M/sec ( +- 0.00% )
3,301,571,922 branches # 794.736 M/sec ( +- 0.00% )
100,029,270 branch-misses # 3.03% of all branches ( +- 0.00% )
4.154441330 seconds time elapsed ( +- 0.01% )
请注意，分支未命中与重复计数大致相同：内循环每次都在最后预测错误。展开以将循环迭代计数保持在大约 22 次以下将使模式足够短，以便 Skylake 的分支预测器在大多数情况下正确预测未采用的条件。分支错误预测是我们没有通过管道每个周期获得 ~4.0 uop 的唯一原因，因此避免分支未命中会将 IPC 从 3.5 提高到 4.0 以上（cmp/jcc 宏融合将 2 条指令放入一个 uop）。
即使您在 L2 缓存带宽（而不是前端）上遇到瓶颈，这些分支未命中也可能会受到伤害。不过，我没有对此进行测试：我的测试只是for()围绕来自 Paul R 的测试工具的函数调用进行了循环，因此 L1D 缓存中的所有内容都很热门。内部循环的 32 次迭代接近这里的最坏情况：足够低以防止频繁的错误预测，但不会低到分支预测可以获取模式并避免它们。
我的版本应该在每次迭代中运行 3 个周期，仅在前端、英特尔 Sandybridge 及更高版本上遇到瓶颈。（Nehalem 将在每个时钟一个负载上出现瓶颈。）
参见http://agner.org/optimize/，还有x86 的 MOV 真的可以“免费”吗？为什么我完全不能重现这个？有关融合域与未融合域 uops 和性能计数器的更多信息。
更新： clang 为您展开它，至少当大小是编译时常量时......奇怪的是，它甚至展开dst[i/2]函数的非内联版本（未知size），但不是较低循环开销的版本。
使用clang++-4.0 -O3 -march=skylake -mno-avx，我的版本（由编译器展开 2）运行： 9.61G 周期，100M 迭代（2.2 秒）。（35.6G uops 发布（融合域），45.0G uops 执行（未融合域），几乎为零的分支未命中。）可能不再是前端的瓶颈，但 AVX 仍然会受到伤害。
Paul R's（也由 2 展开）以 12.29G 周期运行 100M 迭代（2.8 秒）。发出 48.4G uops（融合域），执行 51.4G uops（未融合域）。50.1G 指令，对于 4.08 IPC，可能仍然在前端遇到瓶颈（因为它需要几条movdqa指令在销毁寄存器之前复制它）。AVX 将有助于非破坏性向量指令，即使没有 AVX2 用于更宽的整数向量。
通过仔细编码，您应该能够很好地处理运行时变量的大小。
使用对齐的指针和对齐的加载，因此编译器可以使用pavgb内存操作数，而不是使用单独的未对齐加载指令。这意味着前端的指令和 uops 更少，这是此循环的瓶颈。
这对 Paul 的版本没有帮助，因为只有 for 的第二个操作数pmaddubsw可以来自内存，而且它被视为有符号字节。如果我们使用_mm_maddubs_epi16(_mm_set1_epi8(1), v0);，则 16 位乘法结果将进行符号扩展而不是零扩展。所以1+255会出现 0 而不是 256。
折叠负载需要与 SSE 对齐，但不需要与 AVX 对齐。然而，在 Intel Haswell/Skylake 上，索引寻址模式只能与读取-修改-写入目标寄存器的指令保持微融合。 vpavgb xmm0, xmm0, [rsi+rax*2]在进入核心的乱序部分之前，在 Haswell/Skylake 上未层压到 2 uop，但pavgb xmm1, [rsi+rax*2]可以一直保持微融合，因此它作为单个 uop 发出。除了 Ryzen（即不是 Atom/Silvermont），在主流 x86 CPU 上，前端问题的瓶颈是每个时钟 4 个融合域 uops。将一半的负载折叠到内存操作数中有助于在除 Sandybridge/Ivybridge 之外的所有 Intel CPU 和所有 AMD CPU 上实现这一点。
内联到使用 的测试函数时，gcc 和 clang 会折叠负载alignas(32)，即使您使用_mm_loadu内在函数。他们知道数据是一致的，并利用。
奇怪的事实：在启用 AVX 代码生成的情况下编译 128b 向量化代码 ( -march=native) 实际上会在 Haswell/Skylake 上减慢它的速度，因为即使它们是 的内存操作数vpavgb，它也会使所有 4 个负载作为单独的 uops 发出，并且没有t 任何movdqaAVX 会避免的寄存器复制指令。（通常 AVX 无论如何都会领先，即使对于仍然只使用 128b 向量的手动向量化代码，因为 3 操作数指令的好处是不会破坏它们的输入之一。）在这种情况下，13,53G cycles ( +- 0.05% )或者3094.195773 ms ( +- 0.20% )，从11.92G大约 2.7 秒的周期开始. uops_issued = 48.508G，从41,908. 指令计数和 uops_executed 计数本质上是相同的。
OTOH，实际的 256b AVX2 版本的运行速度略低于两倍。一些展开以减少前端瓶颈肯定会有所帮助。根据@Mysticial 的测试，AVX512 版本在 Skylake-AVX512 Xeon 上的运行速度可能接近 4 倍，但可能会导致 ALU 吞吐量瓶颈，因为当 RS 中有任何 512b uops 等待执行时，SKX 会关闭执行端口 1。（这就解释了为什么pavgb zmm每个时钟有 1 个吞吐量，而pavgb ymm每个时钟有 2 个吞吐量。.）
要对齐两个输入行，请以行间距为 16 倍数的格式存储图像数据，即使实际图像尺寸为奇数。您的存储步幅不必与您的实际图像尺寸相匹配。
如果您只能对齐源或目标（例如，因为您要缩小从源图像中奇数列开始的区域），您可能仍然应该对齐源指针。
英特尔的优化手册建议对齐目标而不是源，如果您不能同时对齐两者，但执行 4 倍于存储的负载可能会改变平衡。
要在开始/结束处处理未对齐，请从开始和结束处做一个潜在重叠的未对齐像素向量。商店可以与其他商店重叠，并且由于 dst 与 src 是分开的，因此您可以重做部分重叠的向量。
在 Paul 的 test 中main()，我只是alignas(32)在每个数组前面添加。
由于您使用编译一个版本-march=native，因此您可以在编译时使用 轻松检测 AVX2 #ifdef __AVX2__。没有简单的方法可以对 128b 和 256b 手动矢量化使用完全相同的代码。所有内在函数都有不同的名称，因此即使没有其他差异，您通常也需要复制所有内容。
（有一些 C++ 包装库用于使用运算符重载和函数重载的内在函数，让您编写一个模板化版本，在不同宽度的向量上使用相同的逻辑。例如 Agner Fog 的 VCL 很好，但除非您的软件是开放的 -源，您不能使用它，因为它是 GPL 许可的，并且您想分发二进制文件。）
要在二进制分发版本中利用 AVX2，您必须进行运行时检测/分发。在这种情况下，您希望分派到循环行的函数的版本，因此您在循环行内没有分派开销。或者让那个版本使用 SSSE3。
| 归档时间： |
|
| 查看次数： |
750 次 |
| 最近记录： |