使用 SSE 最快缩小 8 位灰度图像
Fastest downscaling of 8bit gray image with SSE
我有一个函数可以将 8 位图像缩小两倍。我以前 optimised the rgb32 case with SSE。现在我想对 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));
}
}
备注：
- 我意识到我的函数有轻微的（相差一个）舍入错误，但我愿意接受这一点。
- 为清楚起见，我假设大小是 32 的倍数。
编辑：现在有一个 github repository implementing the answers to this question. The fastest solution was provided by user 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);
}
这是一个使用较少指令的实现。不过我还没有根据你的代码对它进行基准测试，所以它可能不会快得多：
void average2Rows(const uint8_t* src1, const uint8_t* src2, uint8_t* dst, int size)
{
const __m128i vk1 = _mm_set1_epi8(1);
for (int i = 0; i < size - 31; i += 32)
{
__m128i v0 = _mm_loadu_si128((__m128i *)&src1[i]);
__m128i v1 = _mm_loadu_si128((__m128i *)&src1[i + 16]);
__m128i v2 = _mm_loadu_si128((__m128i *)&src2[i]);
__m128i v3 = _mm_loadu_si128((__m128i *)&src2[i + 16]);
__m128i w0 = _mm_maddubs_epi16(v0, vk1); // unpack and horizontal add
__m128i w1 = _mm_maddubs_epi16(v1, vk1);
__m128i w2 = _mm_maddubs_epi16(v2, vk1);
__m128i w3 = _mm_maddubs_epi16(v3, vk1);
w0 = _mm_add_epi16(w0, w2); // vertical add
w1 = _mm_add_epi16(w1, w3);
w0 = _mm_srli_epi16(w0, 2); // divide by 4
w1 = _mm_srli_epi16(w1, 2);
w0 = _mm_packus_epi16(w0, w1); // pack
_mm_storeu_si128((__m128i *)&dst[i / 2], w0);
}
}
测试工具：
#include <stdio.h>
#include <stdlib.h>
#include <tmmintrin.h>
void average2Rows_ref(const uint8_t* row1, const uint8_t* row2, uint8_t* dst, int size)
{
for (int i = 0; i < size - 1; i += 2)
{
dst[i / 2] = (row1[i] + row1[i + 1] + row2[i] + row2[i + 1]) / 4;
}
}
void average2Rows(const uint8_t* src1, const uint8_t* src2, uint8_t* dst, int size)
{
const __m128i vk1 = _mm_set1_epi8(1);
for (int i = 0; i < size - 31; i += 32)
{
__m128i v0 = _mm_loadu_si128((__m128i *)&src1[i]);
__m128i v1 = _mm_loadu_si128((__m128i *)&src1[i + 16]);
__m128i v2 = _mm_loadu_si128((__m128i *)&src2[i]);
__m128i v3 = _mm_loadu_si128((__m128i *)&src2[i + 16]);
__m128i w0 = _mm_maddubs_epi16(v0, vk1); // unpack and horizontal add
__m128i w1 = _mm_maddubs_epi16(v1, vk1);
__m128i w2 = _mm_maddubs_epi16(v2, vk1);
__m128i w3 = _mm_maddubs_epi16(v3, vk1);
w0 = _mm_add_epi16(w0, w2); // vertical add
w1 = _mm_add_epi16(w1, w3);
w0 = _mm_srli_epi16(w0, 2); // divide by 4
w1 = _mm_srli_epi16(w1, 2);
w0 = _mm_packus_epi16(w0, w1); // pack
_mm_storeu_si128((__m128i *)&dst[i / 2], w0);
}
}
int main()
{
const int n = 1024;
uint8_t src1[n];
uint8_t src2[n];
uint8_t dest_ref[n / 2];
uint8_t dest_test[n / 2];
for (int i = 0; i < n; ++i)
{
src1[i] = rand();
src2[i] = rand();
}
for (int i = 0; i < n / 2; ++i)
{
dest_ref[i] = 0xaa;
dest_test[i] = 0x55;
}
average2Rows_ref(src1, src2, dest_ref, n);
average2Rows(src1, src2, dest_test, n);
for (int i = 0; i < n / 2; ++i)
{
if (dest_test[i] != dest_ref[i])
{
printf("%u %u %u %u: ref = %u, test = %u\n", src1[2 * i], src1[2 * i + 1], src2[2 * i], src2[2 * i + 1], dest_ref[i], dest_test[i]);
}
}
return 0;
}
请注意，SIMD 版本的输出与标量参考代码的输出完全匹配（没有 "off by one" 舍入错误）。
如果你愿意接受两次使用 pavgb 的双舍入，你可以通过先用 pavgb 进行垂直平均来比 Paul R 的回答更快，减少一半的数量需要解压缩为 16 位元素的数据。 （并允许一半的负载折叠到 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) 将每个水平对的奇数元素与偶数元素对齐。但是因为没有带 t运行cation 的水平包装，所以在包装之前你必须 _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 微指令，它只需要 1 个输入并产生半角输出。根据 IACA，内存目标形式总共有 4 个未融合域 uops（与 reg,reg + separate store 相同）。我不得不使用 _mm_mask_cvtepi16_storeu_epi8(&dst\[i+0\], -1, l_avg); 来获取它，因为 gcc 和 clang 无法将单独的 _mm_store 折叠到 vpmovwb 的内存目标中。 （没有非屏蔽存储内在函数，因为编译器应该像将 _mm_load 折叠到典型 ALU 指令的内存操作数中那样为您做这件事）。
它可能仅在缩小到 1/4 或 1/8 (cvtepi64_epi8) 时有用，而不仅仅是缩小到一半。或者可能有助于避免需要第二次洗牌来处理 _mm512_packus_epi16 的车道内行为。使用 AVX2，在 [D C] [B A] 上的 _mm256_packus_epi16 之后，您有 [D B | C A]，您可以使用 AVX2 _mm256_permute4x64_epi64 (__m256i a, const int imm8) 修复它以在 64 位块中随机播放。但是对于 AVX512，您需要 vpermq 的向量随机播放控件。 packus + 修复洗牌可能仍然是更好的选择。
一旦你这样做了，循环中就没有多少向量指令了，让编译器更紧凑的 asm 有很多好处。不幸的是，你的循环很难让编译器做好。
（这也有助于 Paul R 的解决方案，因为他从问题中复制了对编译器不友好的循环结构。）
以gcc/clang可以更好地优化的方式使用循环计数器，并使用避免每次通过循环重新进行符号扩展的类型。
在您当前的循环中，gcc/clang 实际上对 i/2 进行了算术右移，而不是递增 16（而不是 32）并为负载使用缩放索引寻址模式。他们似乎没有意识到 i 总是偶数。
(full code + asm on Matt Godbolt's compiler explorer):
.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 的指针增量，以及存储的 counter/index。 （完全脑残的行为，特别是因为他们仍然用 -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 次重复，这相差 ~2.7 秒，而 ~3.3 秒。
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 微指令的唯一原因，因此避免分支未命中会将 IPC 从 3.5 提高到 4.0 以上（cmp/jcc 宏融合将 2 条指令放入一个微指令中）。
即使您在 L2 缓存带宽（而不是前端）上遇到瓶颈，这些分支未命中也可能会受到伤害。不过，我没有对此进行测试：我的测试只是围绕来自 Paul R 的测试工具的函数调用包装了一个 for() 循环，因此 L1D 缓存中的所有内容都是热的。内部循环的 32 次迭代接近这里的最坏情况：对于频繁的错误预测来说足够低，但又不会低到分支预测可以选择模式并避免它们。
我的版本应该 运行 每次迭代 3 个周期，仅在前端、Intel Sandybridge 和更高版本上出现瓶颈。 （Nehalem 将在每个时钟一个负载上出现瓶颈。）
有关融合域与非融合域微指令和性能计数器的更多信息，请参阅 http://agner.org/optimize/, and also 。
更新：clang 为您展开它，至少当大小是编译时常量时......奇怪的是，它甚至展开 dst[i/2] 函数（未知 size），但不是低循环开销版本。
使用 clang++-4.0 -O3 -march=skylake -mno-avx，我的版本（由编译器展开 2）运行s：100M 迭代（2.2s）的 9.61G 周期。 （发布 35.6G uops（融合域），执行 45.0G uops（非融合域），分支未命中率接近零。）可能不再是前端瓶颈，但 AVX 仍然会受到伤害。
Paul R 的（也展开了 2）运行s 在 12.29G 周期中用于 100M 迭代（2.8s）。发出 48.4G uops（融合域），执行 51.4G uops（非融合域）。 50.1G 指令，对于 4.08 IPC，可能仍然是前端的瓶颈（因为它需要一对 movdqa 指令来在销毁寄存器之前复制寄存器）。 AVX 将有助于非破坏性向量指令，即使没有用于更宽整数向量的 AVX2。
通过仔细编码，您应该能够很好地处理 运行 时变大小。
使用对齐指针和对齐加载，因此编译器可以使用 pavgb 和内存操作数 而不是使用单独的未对齐加载指令。这意味着前端的指令和微指令更少，这是此循环的瓶颈。
这对 Paul 的版本没有帮助，因为只有 pmaddubsw 的第二个操作数可以来自内存，也就是被视为有符号字节的那个。如果我们使用 _mm_maddubs_epi16(_mm_set1_epi8(1), v0);，16 位乘法结果将被符号扩展而不是零扩展。所以 1+255 会变成 0 而不是 256。
折叠负载需要与 SSE 对齐，但不需要与 AVX 对齐。但是，on Intel Haswell/Skylake, indexed addressing modes can only stay micro-fused with instructions which read-modify-write their destination register. vpavgb xmm0, xmm0, [rsi+rax*2] 在 Haswell/Skylake 上未层压到 2 微指令，然后才进入核心的乱序部分，但 pavgb xmm1, [rsi+rax*2] 可以保持微一路融合，因此它作为单个 uop 发布。前端问题瓶颈是主流 x86 CPU 上每时钟 4 个融合域微指令，除了 Ryzen（即不是 Atom/Silvermont）。将一半的负载折叠到内存操作数中有助于除 Sandybridge/Ivybridge 之外的所有 Intel CPU 和所有 AMD CPU。
gcc 和 clang 在内联到使用 alignas(32) 的测试函数时会折叠负载，即使您使用 _mm_loadu 内部函数也是如此。他们知道数据是一致的，并加以利用。
奇怪的事实：在启用 AVX 代码生成 (-march=native) 的情况下编译 128b 向量化代码实际上会在 Haswell/Skylake 上减慢它的速度，因为它会使所有 4 个加载问题甚至作为单独的 uops当它们是 vpavgb 的内存操作数时，AVX 不会避免任何 movdqa 寄存器复制指令。 （通常 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 Xeons 上的速度可能 运行 接近 4 倍，但可能会成为 ALU 吞吐量的瓶颈，因为当 RS 中有任何 512b 微指令等待执行时，SKX 会关闭执行端口 1 . （这解释了为什么 pavgb zmm has 1 per clock throughput while pavgb ymm is 2 per clock.。）
要使两个输入行对齐，以行跨度为 16 的倍数的格式存储图像数据，即使实际图像尺寸是奇数也是如此。您的存储步幅不必与您的实际图像尺寸相匹配。
如果您只能对齐源指针或目标指针（例如，因为您正在缩小从源图像中的奇数列开始的区域），您可能仍应对齐源指针。
Intel 的优化手册建议对齐目的地而不是源，如果您不能对齐两者，但是做 4 倍于存储的负载可能会改变平衡。
要处理 start/end 处的未对齐，从开始和结束处执行可能重叠的未对齐像素向量。 store 与其他 store 重叠是可以的，并且由于 dst 与 src 是分开的，因此您可以重做一个部分重叠的向量。
在保罗的测试main()中，我只是在每个数组前面添加了alignas(32)。
AVX2:
因为你 ，你可以在编译时使用 #ifdef __AVX2__ 轻松检测 AVX2。没有简单的方法可以为 128b 和 256b 手动矢量化使用完全相同的代码。所有内在函数都有不同的名称，因此即使没有其他差异，您通常也需要复制所有内容。
（有一些使用运算符重载和函数重载的内在函数的 C++ 包装器库，让您编写在不同宽度的向量上使用相同逻辑的模板化版本。例如，Agner Fog 的 VCL 很好，但除非你的软件是开源的，你不能使用它，因为它是 GPL 许可的，你想分发一个二进制文件。）
要在您的二进制分发版本中利用 AVX2，您必须 运行time detection/dispatching。在这种情况下，您希望分派到循环遍历行的函数版本，这样您就不会在遍历行的循环中产生分派开销。或者让那个版本使用 SSSE3。
我有一个函数可以将 8 位图像缩小两倍。我以前 optimised the rgb32 case with SSE。现在我想对 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));
}
}
备注：
- 我意识到我的函数有轻微的（相差一个）舍入错误，但我愿意接受这一点。
- 为清楚起见，我假设大小是 32 的倍数。
编辑：现在有一个 github repository implementing the answers to this question. The fastest solution was provided by user 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);
}
这是一个使用较少指令的实现。不过我还没有根据你的代码对它进行基准测试，所以它可能不会快得多：
void average2Rows(const uint8_t* src1, const uint8_t* src2, uint8_t* dst, int size)
{
const __m128i vk1 = _mm_set1_epi8(1);
for (int i = 0; i < size - 31; i += 32)
{
__m128i v0 = _mm_loadu_si128((__m128i *)&src1[i]);
__m128i v1 = _mm_loadu_si128((__m128i *)&src1[i + 16]);
__m128i v2 = _mm_loadu_si128((__m128i *)&src2[i]);
__m128i v3 = _mm_loadu_si128((__m128i *)&src2[i + 16]);
__m128i w0 = _mm_maddubs_epi16(v0, vk1); // unpack and horizontal add
__m128i w1 = _mm_maddubs_epi16(v1, vk1);
__m128i w2 = _mm_maddubs_epi16(v2, vk1);
__m128i w3 = _mm_maddubs_epi16(v3, vk1);
w0 = _mm_add_epi16(w0, w2); // vertical add
w1 = _mm_add_epi16(w1, w3);
w0 = _mm_srli_epi16(w0, 2); // divide by 4
w1 = _mm_srli_epi16(w1, 2);
w0 = _mm_packus_epi16(w0, w1); // pack
_mm_storeu_si128((__m128i *)&dst[i / 2], w0);
}
}
测试工具：
#include <stdio.h>
#include <stdlib.h>
#include <tmmintrin.h>
void average2Rows_ref(const uint8_t* row1, const uint8_t* row2, uint8_t* dst, int size)
{
for (int i = 0; i < size - 1; i += 2)
{
dst[i / 2] = (row1[i] + row1[i + 1] + row2[i] + row2[i + 1]) / 4;
}
}
void average2Rows(const uint8_t* src1, const uint8_t* src2, uint8_t* dst, int size)
{
const __m128i vk1 = _mm_set1_epi8(1);
for (int i = 0; i < size - 31; i += 32)
{
__m128i v0 = _mm_loadu_si128((__m128i *)&src1[i]);
__m128i v1 = _mm_loadu_si128((__m128i *)&src1[i + 16]);
__m128i v2 = _mm_loadu_si128((__m128i *)&src2[i]);
__m128i v3 = _mm_loadu_si128((__m128i *)&src2[i + 16]);
__m128i w0 = _mm_maddubs_epi16(v0, vk1); // unpack and horizontal add
__m128i w1 = _mm_maddubs_epi16(v1, vk1);
__m128i w2 = _mm_maddubs_epi16(v2, vk1);
__m128i w3 = _mm_maddubs_epi16(v3, vk1);
w0 = _mm_add_epi16(w0, w2); // vertical add
w1 = _mm_add_epi16(w1, w3);
w0 = _mm_srli_epi16(w0, 2); // divide by 4
w1 = _mm_srli_epi16(w1, 2);
w0 = _mm_packus_epi16(w0, w1); // pack
_mm_storeu_si128((__m128i *)&dst[i / 2], w0);
}
}
int main()
{
const int n = 1024;
uint8_t src1[n];
uint8_t src2[n];
uint8_t dest_ref[n / 2];
uint8_t dest_test[n / 2];
for (int i = 0; i < n; ++i)
{
src1[i] = rand();
src2[i] = rand();
}
for (int i = 0; i < n / 2; ++i)
{
dest_ref[i] = 0xaa;
dest_test[i] = 0x55;
}
average2Rows_ref(src1, src2, dest_ref, n);
average2Rows(src1, src2, dest_test, n);
for (int i = 0; i < n / 2; ++i)
{
if (dest_test[i] != dest_ref[i])
{
printf("%u %u %u %u: ref = %u, test = %u\n", src1[2 * i], src1[2 * i + 1], src2[2 * i], src2[2 * i + 1], dest_ref[i], dest_test[i]);
}
}
return 0;
}
请注意，SIMD 版本的输出与标量参考代码的输出完全匹配（没有 "off by one" 舍入错误）。
如果你愿意接受两次使用 pavgb 的双舍入，你可以通过先用 pavgb 进行垂直平均来比 Paul R 的回答更快，减少一半的数量需要解压缩为 16 位元素的数据。 （并允许一半的负载折叠到 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) 将每个水平对的奇数元素与偶数元素对齐。但是因为没有带 t运行cation 的水平包装，所以在包装之前你必须 _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 微指令，它只需要 1 个输入并产生半角输出。根据 IACA，内存目标形式总共有 4 个未融合域 uops（与 reg,reg + separate store 相同）。我不得不使用 _mm_mask_cvtepi16_storeu_epi8(&dst\[i+0\], -1, l_avg); 来获取它，因为 gcc 和 clang 无法将单独的 _mm_store 折叠到 vpmovwb 的内存目标中。 （没有非屏蔽存储内在函数，因为编译器应该像将 _mm_load 折叠到典型 ALU 指令的内存操作数中那样为您做这件事）。
它可能仅在缩小到 1/4 或 1/8 (cvtepi64_epi8) 时有用，而不仅仅是缩小到一半。或者可能有助于避免需要第二次洗牌来处理 _mm512_packus_epi16 的车道内行为。使用 AVX2，在 [D C] [B A] 上的 _mm256_packus_epi16 之后，您有 [D B | C A]，您可以使用 AVX2 _mm256_permute4x64_epi64 (__m256i a, const int imm8) 修复它以在 64 位块中随机播放。但是对于 AVX512，您需要 vpermq 的向量随机播放控件。 packus + 修复洗牌可能仍然是更好的选择。
一旦你这样做了，循环中就没有多少向量指令了，让编译器更紧凑的 asm 有很多好处。不幸的是，你的循环很难让编译器做好。 （这也有助于 Paul R 的解决方案，因为他从问题中复制了对编译器不友好的循环结构。）
以gcc/clang可以更好地优化的方式使用循环计数器，并使用避免每次通过循环重新进行符号扩展的类型。
在您当前的循环中，gcc/clang 实际上对 i/2 进行了算术右移，而不是递增 16（而不是 32）并为负载使用缩放索引寻址模式。他们似乎没有意识到 i 总是偶数。
(full code + asm on Matt Godbolt's compiler explorer):
.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 的指针增量，以及存储的 counter/index。 （完全脑残的行为，特别是因为他们仍然用 -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 次重复，这相差 ~2.7 秒，而 ~3.3 秒。
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 微指令的唯一原因，因此避免分支未命中会将 IPC 从 3.5 提高到 4.0 以上（cmp/jcc 宏融合将 2 条指令放入一个微指令中）。
即使您在 L2 缓存带宽（而不是前端）上遇到瓶颈，这些分支未命中也可能会受到伤害。不过，我没有对此进行测试：我的测试只是围绕来自 Paul R 的测试工具的函数调用包装了一个 for() 循环，因此 L1D 缓存中的所有内容都是热的。内部循环的 32 次迭代接近这里的最坏情况：对于频繁的错误预测来说足够低，但又不会低到分支预测可以选择模式并避免它们。
我的版本应该 运行 每次迭代 3 个周期，仅在前端、Intel Sandybridge 和更高版本上出现瓶颈。 （Nehalem 将在每个时钟一个负载上出现瓶颈。）
有关融合域与非融合域微指令和性能计数器的更多信息，请参阅 http://agner.org/optimize/, and also
更新：clang 为您展开它，至少当大小是编译时常量时......奇怪的是，它甚至展开 dst[i/2] 函数（未知 size），但不是低循环开销版本。
使用 clang++-4.0 -O3 -march=skylake -mno-avx，我的版本（由编译器展开 2）运行s：100M 迭代（2.2s）的 9.61G 周期。 （发布 35.6G uops（融合域），执行 45.0G uops（非融合域），分支未命中率接近零。）可能不再是前端瓶颈，但 AVX 仍然会受到伤害。
Paul R 的（也展开了 2）运行s 在 12.29G 周期中用于 100M 迭代（2.8s）。发出 48.4G uops（融合域），执行 51.4G uops（非融合域）。 50.1G 指令，对于 4.08 IPC，可能仍然是前端的瓶颈（因为它需要一对 movdqa 指令来在销毁寄存器之前复制寄存器）。 AVX 将有助于非破坏性向量指令，即使没有用于更宽整数向量的 AVX2。
通过仔细编码，您应该能够很好地处理 运行 时变大小。
使用对齐指针和对齐加载，因此编译器可以使用 pavgb 和内存操作数 而不是使用单独的未对齐加载指令。这意味着前端的指令和微指令更少，这是此循环的瓶颈。
这对 Paul 的版本没有帮助，因为只有 pmaddubsw 的第二个操作数可以来自内存，也就是被视为有符号字节的那个。如果我们使用 _mm_maddubs_epi16(_mm_set1_epi8(1), v0);，16 位乘法结果将被符号扩展而不是零扩展。所以 1+255 会变成 0 而不是 256。
折叠负载需要与 SSE 对齐，但不需要与 AVX 对齐。但是，on Intel Haswell/Skylake, indexed addressing modes can only stay micro-fused with instructions which read-modify-write their destination register. vpavgb xmm0, xmm0, [rsi+rax*2] 在 Haswell/Skylake 上未层压到 2 微指令，然后才进入核心的乱序部分，但 pavgb xmm1, [rsi+rax*2] 可以保持微一路融合，因此它作为单个 uop 发布。前端问题瓶颈是主流 x86 CPU 上每时钟 4 个融合域微指令，除了 Ryzen（即不是 Atom/Silvermont）。将一半的负载折叠到内存操作数中有助于除 Sandybridge/Ivybridge 之外的所有 Intel CPU 和所有 AMD CPU。
gcc 和 clang 在内联到使用 alignas(32) 的测试函数时会折叠负载，即使您使用 _mm_loadu 内部函数也是如此。他们知道数据是一致的，并加以利用。
奇怪的事实：在启用 AVX 代码生成 (-march=native) 的情况下编译 128b 向量化代码实际上会在 Haswell/Skylake 上减慢它的速度，因为它会使所有 4 个加载问题甚至作为单独的 uops当它们是 vpavgb 的内存操作数时，AVX 不会避免任何 movdqa 寄存器复制指令。 （通常 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 Xeons 上的速度可能 运行 接近 4 倍，但可能会成为 ALU 吞吐量的瓶颈，因为当 RS 中有任何 512b 微指令等待执行时，SKX 会关闭执行端口 1 . （这解释了为什么 pavgb zmm has 1 per clock throughput while pavgb ymm is 2 per clock.。）
要使两个输入行对齐，以行跨度为 16 的倍数的格式存储图像数据，即使实际图像尺寸是奇数也是如此。您的存储步幅不必与您的实际图像尺寸相匹配。
如果您只能对齐源指针或目标指针（例如，因为您正在缩小从源图像中的奇数列开始的区域），您可能仍应对齐源指针。
Intel 的优化手册建议对齐目的地而不是源，如果您不能对齐两者，但是做 4 倍于存储的负载可能会改变平衡。
要处理 start/end 处的未对齐，从开始和结束处执行可能重叠的未对齐像素向量。 store 与其他 store 重叠是可以的，并且由于 dst 与 src 是分开的，因此您可以重做一个部分重叠的向量。
在保罗的测试main()中，我只是在每个数组前面添加了alignas(32)。
AVX2:
因为你 #ifdef __AVX2__ 轻松检测 AVX2。没有简单的方法可以为 128b 和 256b 手动矢量化使用完全相同的代码。所有内在函数都有不同的名称，因此即使没有其他差异，您通常也需要复制所有内容。
（有一些使用运算符重载和函数重载的内在函数的 C++ 包装器库，让您编写在不同宽度的向量上使用相同逻辑的模板化版本。例如，Agner Fog 的 VCL 很好，但除非你的软件是开源的，你不能使用它，因为它是 GPL 许可的，你想分发一个二进制文件。）
要在您的二进制分发版本中利用 AVX2，您必须 运行time detection/dispatching。在这种情况下，您希望分派到循环遍历行的函数版本，这样您就不会在遍历行的循环中产生分派开销。或者让那个版本使用 SSSE3。