optimization - 使用 AVX2 指令选择性地异或列表的元素
问题描述
我想用 AVX2 指令加速以下操作，但我找不到这样做的方法。
我得到了一个大uint64_t data[100000]的 uint64_t 数组和一个unsigned char indices[100000]字节数组。我想输出一个数组uint64_t Out[256]，其中第 i 个值是所有data[j]这样的 xor index[j]=i。
我想要的一个简单的实现是这样的：
uint64_t Out[256] = {0}; // initialize output array
for (i = 0; i < 100000 ; i++) {
Out[Indices[i]] ^= data[i];
}
我们可以使用 AVX2 指令更有效地实现这一点吗？
编辑：这就是我的代码现在的样子
uint64_t Out[256][4] = {0}; // initialize output array
for (i = 0; i < 100000 ; i+=4) {
Out[Indices[i ]][0] ^= data[i];
Out[Indices[i+1]][1] ^= data[i+1];
Out[Indices[i+2]][2] ^= data[i+2];
Out[Indices[i+3]][3] ^= data[i+3];
}
解决方案
基于对 Haswell/Skylake 的静态分析，我提出了一个版本i，当由 gcc 编译时，每 4 个值运行约 5 个周期，而不是 8 个周期。大尺寸的平均值，不包括组合 的多个副本的时间Out[]，并假设索引的随机分布不会导致任何存储/重新加载 dep 链运行足够长的时间而不重要。
IDK 如果您关心 Ryzen 或 Excavator（其他 2 个主流 AVX2 微架构）。
我没有手工进行仔细分析，但IACA对 HSW/SKL 是错误的，并且认为某些指令实际上不会微熔（在带有 perf 计数器的 i7-6700k 上测试），所以它认为前端瓶颈比实际情况更严重。例如movhps加载+合并微熔断器，但 IACA 认为它甚至不适用于简单的寻址模式。
我们应该可以忽略任何缓存未命中，因为uint64_t Out[4][256]只有 8kiB。因此，在最新的 CPU 上，我们的缓存占用空间仅为 L1d 大小的 1/4，即使在两个逻辑线程之间共享 L1d 的超线程也应该基本没问题。循环data[]并Indices[]应该很好地预取，并且希望不会驱逐Out[]太多。因此，静态分析很可能会比较准确，而且它比仔细的微基准测试更快，更重要的是可以准确地告诉您瓶颈是什么。
但当然，我们严重依赖乱序执行和不完美的调度或其他意想不到的瓶颈很容易发生。不过，如果我没有得到报酬，我并不想真正进行微基准测试。
我们可以使用 AVX2 指令更有效地实现这一点吗？
这基本上是一个直方图问题。使用多个表并在最后组合的通常的直方图优化适用。SIMD XOR 对于最后的组合很有用（只要您使用Out[4][256]，而不是Out[256][4]。后者还通过要求缩放8*4而不是使索引变慢8（这可以通过缩放索引寻址中的单个 LEA 来完成）模式））。
但与正常的直方图不同，您是在内存中对一些数据进行异或运算，而不是添加一个常数 1。因此1，代码必须加载data[i]到寄存器中，而不是立即数，作为xor. （或加载，然后xor reg, data[i]/存储）。这比直方图还要多的总内存操作。
movq我们从“手动”聚集/分散到 SIMD 向量（使用/ loads/stores）中脱颖而出movhps，允许我们使用 SIMD 进行data[i]加载和 XOR。这减少了加载操作的总数，从而减少了加载端口压力，而不会花费额外的前端带宽。
手动收集到 256 位向量可能不值得额外的改组（额外的 vinserti128 / vextracti128 只是为了我们可以将 2 个 memory-source 组合vpxor成一个 256 位的）。128 位向量应该很好。前端吞吐量也是一个主要问题，因为（在 Intel SnB 系列 CPU 上）您希望避免存储的索引寻址模式。gcc 使用lea指令来计算寄存器中的地址，而不是使用索引加载/存储。clang / LLVM-march=skylake决定不这样做，在这种情况下这是一个糟糕的决定，因为端口 2 / 端口 3 上的循环瓶颈，并且花费额外的 ALU 微指令来允许存储地址微指令使用端口 7 是一个胜利。但是，如果您在 p23 上没有遇到瓶颈，那么花费额外的微指令来避免索引存储是不好的。（和在可以保持微融合的情况下，绝对不仅仅是为了避免索引负载；愚蠢的GCC）。也许 gcc 和 LLVM 的寻址模式成本模型不是很准确，或者它们没有对管道进行足够详细的建模，无法确定前端与特定端口的循环瓶颈何时出现。
寻址模式的选择和其他 asm 代码生成选择对于在 SnB 系列上实现最佳性能至关重要。但是用 C 语言编写使您无法控制它；您主要受编译器的支配，除非您可以调整源代码以使其做出不同的选择。例如 gcc 与 clang 在这里有很大的不同。
在 SnB 系列上，movhps负载需要端口 5 用于 shuffle/blend（尽管它确实将微融合到一个 uop 中），但movhpsstore 是没有 ALU uop 的纯 store。所以它是收支平衡的，让我们对两个数据元素使用一个 SIMD 加载/异或。
使用 AVX，ALU 微指令允许未对齐的内存源操作数，因此我们不需要为data[]. 但英特尔 HSW/SKL 可以保持索引寻址模式微融合pxor但不能vpxor。所以在不启用 AVX 的情况下编译会更好，允许编译器使用索引寻址模式而不是递增单独的指针。（或者，如果编译器不知道这一点并且无论如何都使用索引寻址模式，则使其更快。） TL：DR：可能最好要求 16 字节对齐data[]并在禁用 AVX 的情况下编译该函数，以获得更好的宏融合。（但是我们错过了在最后组合切片的 256 位 SIMD Out，除非我们把它放在用 AVX 或 AVX2 编译的不同函数中）
避免未对齐的负载也将避免任何缓存行拆分，这不会花费额外的微指令，但我们可能接近 L1d 吞吐量限制的瓶颈，而不仅仅是加载/存储执行单元的吞吐量限制。
我还查看了一次加载 4 个索引并使用 ALU 指令解包。例如与memcpyinto struct { uint8_t idx[4]; } idx;。但是 gcc 会生成多个用于解包的浪费指令。太糟糕了 x86 没有像 ARMubfx或特别是PowerPC这样的出色位域指令rlwinm（这可能会使结果免费左移，所以如果 x86 有这个，静态Out可能会在非 PIC 代码中使用 base+disp32 寻址模式。 )
如果我们使用标量 XOR，则使用 AL/AH 中的 shift / movzx 解包 dword 是一个胜利，但看起来不是当我们使用 SIMDdata[]并将前端吞吐量用于lea指令以允许存储地址 uops 时在端口 7 上运行。这使我们成为前端瓶颈而不是端口 2/3 的瓶颈，因此movzx根据静态分析，使用 4 倍内存负载看起来最好。如果您花时间手动编辑 asm，那么两种方式都值得进行基准测试。（gcc 生成的带有额外微指令的 asm 很糟糕，包括movzx在右移 24 位后完全冗余，高位已经为零。）
编码
（在 Godbolt 编译器资源管理器中查看它，以及标量版本）：
#include <immintrin.h>
#include <stdint.h>
#include <string.h>
#include <stdalign.h>
#ifdef IACA_MARKS
#include "/opt/iaca-3.0/iacaMarks.h"
#else
#define IACA_START
#define IACA_END
#endif
void hist_gatherscatter(unsigned idx0, unsigned idx1,
uint64_t Out0[256], uint64_t Out1[256],
__m128i vdata) {
// gather load from Out[0][?] and Out[1][?] with movq / movhps
__m128i hist = _mm_loadl_epi64((__m128i*)&Out0[idx0]);
hist = _mm_castps_si128( // movhps into the high half
_mm_loadh_pi(_mm_castsi128_ps(hist), (__m64*)&Out1[idx1]));
// xorps could bottleneck on port5.
// Actually probably not, using __m128 the whole time would be simpler and maybe not confuse clang
hist = _mm_xor_si128(hist, vdata);
// scatter store with movq / movhps
_mm_storel_epi64((__m128i*)&Out0[idx0], hist);
_mm_storeh_pi((__m64*)&Out1[idx1], _mm_castsi128_ps(hist));
}
void ext(uint64_t*);
void xor_histo_avx(uint8_t *Indices, const uint64_t *data, size_t len)
{
alignas(32) uint64_t Out[4][256] = {{0}};
// optional: peel the first iteration and optimize away loading the old known-zero values from Out[0..3][Indices[0..3]].
if (len<3) // not shown: cleanup for last up-to-3 elements.
return;
for (size_t i = 0 ; i<len ; i+=4) {
IACA_START
// attempt to hand-hold compiler into a dword load + shifts to extract indices
// to reduce load-port pressure
struct { uint8_t idx[4]; } idx;
#if 0
memcpy(&idx, Indices+i, sizeof(idx)); // safe with strict-aliasing and possibly-unaligned
//gcc makes stupid asm for this, same as for memcpy into a struct,
// using a dword load into EAX (good),
// then AL/AH for the first 2 (good)
// but then redundant mov and movzx instructions for the high 2
// clang turns it into 4 loads
/*
//Attempt to hand-hold gcc into less-stupid asm
//doesn't work: same asm as the struct
uint32_t tmp;
memcpy(&tmp, Indices+i, sizeof(tmp)); // mov eax,[mem]
idx.idx[0] = tmp; //movzx reg, AL
idx.idx[1] = tmp>>8; //movzx reg, AH
tmp >>= 16; //shr eax, 16
idx.idx[2] = tmp; //movzx reg, AL
idx.idx[3] = tmp>>8; //movzx reg, AH
*/
#else
// compiles to separate loads with gcc and clang
idx.idx[0] = Indices[i+0];
idx.idx[1] = Indices[i+1];
idx.idx[2] = Indices[i+2];
idx.idx[3] = Indices[i+3];
#endif
__m128i vd = _mm_load_si128((const __m128i*)&data[i]);
hist_gatherscatter(idx.idx[0], idx.idx[1], Out[0], Out[1], vd);
vd = _mm_load_si128((const __m128i*)&data[i+2]);
hist_gatherscatter(idx.idx[2], idx.idx[3], Out[2], Out[3], vd);
}
IACA_END
// hand-hold compilers into a pointer-increment loop
// to avoid indexed addressing modes. (4/5 speedup on HSW/SKL if all the stores use port7)
__m256i *outp = (__m256i*)&Out[0];
__m256i *endp = (__m256i*)&Out[3][256];
for (; outp < endp ; outp++) {
outp[0] ^= outp[256/4*1];
outp[0] ^= outp[256/4*2];
outp[0] ^= outp[256/4*3];
}
// This part compiles horribly with -mno-avx, but does compile
// because I used GNU C native vector operators on __m256i instead of intrinsics.
/*
for (int i=0 ; i<256 ; i+=4) {
// use loadu / storeu if Out isn't aligned
__m256i out0 = _mm256_load_si256(&Out[0][i]);
__m256i out1 = _mm256_load_si256(&Out[1][i]);
__m256i out2 = _mm256_load_si256(&Out[2][i]);
__m256i out3 = _mm256_load_si256(&Out[3][i]);
out0 = _mm256_xor_si256(out0, out1);
out0 = _mm256_xor_si256(out0, out2);
out0 = _mm256_xor_si256(out0, out3);
_mm256_store_si256(&Out[0][i], out0);
}
*/
//ext(Out[0]); // prevent optimizing away the work
asm("" :: "r"(Out) : "memory");
}
用 gcc7.3 编译，-std=gnu11 -DIACA_MARKS -O3 -march=skylake -mno-avx用 IACA-3.0 分析：
$ /opt/iaca-3.0/iaca xor-histo.iaca.o Intel(R) Architecture Code Analyzer Version - v3.0-28-g1ba2cbb build date: 2017-10-23;16:42:45
Analyzed File - xor-histo.iaca.o
Binary Format - 64Bit
Architecture - SKL
Analysis Type - Throughput
Throughput Analysis Report
--------------------------
Block Throughput: 5.79 Cycles Throughput Bottleneck: FrontEnd
Loop Count: 22 (this is fused-domain uops. It's actually 20, so a 5 cycle front-end bottleneck)
Port Binding In Cycles Per Iteration:
--------------------------------------------------------------------------------------------------
| Port | 0 - DV | 1 | 2 - D | 3 - D | 4 | 5 | 6 | 7 |
--------------------------------------------------------------------------------------------------
| Cycles | 2.0 0.0 | 3.0 | 5.5 5.1 | 5.5 4.9 | 4.0 | 3.0 | 2.0 | 3.0 |
--------------------------------------------------------------------------------------------------
DV - Divider pipe (on port 0)
D - Data fetch pipe (on ports 2 and 3)
F - Macro Fusion with the previous instruction occurred
* - instruction micro-ops not bound to a port
^ - Micro Fusion occurred
# - ESP Tracking sync uop was issued
@ - SSE instruction followed an AVX256/AVX512 instruction, dozens of cycles penalty is expected
X - instruction not supported, was not accounted in Analysis
| Num Of | Ports pressure in cycles | |
| Uops | 0 - DV | 1 | 2 - D | 3 - D | 4 | 5 | 6 | 7 |
-----------------------------------------------------------------------------------------
| 1 | | | 0.5 0.5 | 0.5 0.5 | | | | | movzx r8d, byte ptr [rdi]
| 1 | | | 0.5 0.5 | 0.5 0.5 | | | | | movzx edx, byte ptr [rdi+0x2]
| 1 | | | | | | | 1.0 | | add rdi, 0x4
| 1 | | | | | | | 1.0 | | add rsi, 0x20
| 1 | | | 0.5 0.5 | 0.5 0.5 | | | | | movzx eax, byte ptr [rdi-0x1]
| 1 | | 1.0 | | | | | | | lea r12, ptr [rcx+r8*8]
| 1 | | | 0.5 0.5 | 0.5 0.5 | | | | | movzx r8d, byte ptr [rdi-0x3]
| 1 | | 1.0 | | | | | | | lea rdx, ptr [r10+rdx*8]
| 1 | | | 0.5 0.5 | 0.5 0.5 | | | | | movq xmm0, qword ptr [r12]
| 1 | | | | | | 1.0 | | | lea rax, ptr [r9+rax*8]
| 1 | | 1.0 | | | | | | | lea r8, ptr [r11+r8*8]
| 2 | | | 0.5 0.5 | 0.5 0.5 | | 1.0 | | | movhps xmm0, qword ptr [r8] # Wrong, 1 micro-fused uop on SKL
| 2^ | 1.0 | | 0.5 0.5 | 0.5 0.5 | | | | | pxor xmm0, xmmword ptr [rsi-0x20]
| 2^ | | | 0.5 | 0.5 | 1.0 | | | | movq qword ptr [r12], xmm0 # can run on port 7, IDK why IACA chooses not to model it there
| 2^ | | | | | 1.0 | | | 1.0 | movhps qword ptr [r8], xmm0
| 1 | | | 0.5 0.5 | 0.5 0.5 | | | | | movq xmm0, qword ptr [rdx]
| 2 | | | 0.5 0.5 | 0.5 0.5 | | 1.0 | | | movhps xmm0, qword ptr [rax] # Wrong, 1 micro-fused uop on SKL
| 2^ | 1.0 | | 0.5 0.5 | 0.5 0.5 | | | | | pxor xmm0, xmmword ptr [rsi-0x10]
| 2^ | | | | | 1.0 | | | 1.0 | movq qword ptr [rdx], xmm0
| 2^ | | | | | 1.0 | | | 1.0 | movhps qword ptr [rax], xmm0
| 1* | | | | | | | | | cmp rbx, rdi
| 0*F | | | | | | | | | jnz 0xffffffffffffffa0
Total Num Of Uops: 29 (This is unfused-domain, and a weird thing to total up).
Godbolt 上的 gcc8.1 对 使用缩放索引寻址模式pxor，对 Indices 和 使用相同的计数器data[]，这样可以节省add.
clang 不使用 LEA，并且i每 7 个周期 4 秒的瓶颈，因为没有一个存储 uop 可以在端口 7 上运行。
标量版本（仍然使用 4 片Out[4][256]）：
$ iaca.sh -mark 2 xor-histo.iaca.o
Intel(R) Architecture Code Analyzer Version - 2.3 build:246dfea (Thu, 6 Jul 2017 13:38:05 +0300)
Analyzed File - xor-histo.iaca.o
Binary Format - 64Bit
Architecture - SKL
Analysis Type - Throughput
*******************************************************************
Intel(R) Architecture Code Analyzer Mark Number 2
*******************************************************************
Throughput Analysis Report
--------------------------
Block Throughput: 7.24 Cycles Throughput Bottleneck: FrontEnd
Port Binding In Cycles Per Iteration:
---------------------------------------------------------------------------------------
| Port | 0 - DV | 1 | 2 - D | 3 - D | 4 | 5 | 6 | 7 |
---------------------------------------------------------------------------------------
| Cycles | 3.0 0.0 | 3.0 | 6.2 4.5 | 6.8 4.5 | 4.0 | 3.0 | 3.0 | 0.0 |
---------------------------------------------------------------------------------------
N - port number or number of cycles resource conflict caused delay, DV - Divider pipe (on port 0)
D - Data fetch pipe (on ports 2 and 3), CP - on a critical path
F - Macro Fusion with the previous instruction occurred
* - instruction micro-ops not bound to a port
^ - Micro Fusion happened
# - ESP Tracking sync uop was issued
@ - SSE instruction followed an AVX256/AVX512 instruction, dozens of cycles penalty is expected
X - instruction not supported, was not accounted in Analysis
| Num Of | Ports pressure in cycles | |
| Uops | 0 - DV | 1 | 2 - D | 3 - D | 4 | 5 | 6 | 7 | |
---------------------------------------------------------------------------------
| 1 | | | 0.5 0.5 | 0.5 0.5 | | | | | | mov eax, dword ptr [rdi]
| 1 | 0.4 | 0.5 | | | | 0.1 | | | | add rdi, 0x4
| 1 | | 0.7 | | | | 0.3 | | | | add rsi, 0x20
| 1* | | | | | | | | | | movzx r9d, al
| 1 | | | 0.5 0.5 | 0.5 0.5 | | | | | | mov rdx, qword ptr [rbp+r9*8-0x2040]
| 2^ | | 0.3 | 0.5 0.5 | 0.5 0.5 | | 0.3 | 0.4 | | | xor rdx, qword ptr [rsi-0x20]
| 2 | | | 0.5 | 0.5 | 1.0 | | | | | mov qword ptr [rbp+r9*8-0x2040], rdx # wrong, HSW/SKL can keep indexed stores fused
| 1* | | | | | | | | | | movzx edx, ah
| 1 | | | | | | 0.4 | 0.6 | | | add rdx, 0x100
| 1 | | | 0.5 0.5 | 0.5 0.5 | | | | | | mov r9, qword ptr [rbp+rdx*8-0x2040]
| 2^ | 0.6 | 0.2 | 0.5 0.5 | 0.5 0.5 | | 0.2 | 0.1 | | | xor r9, qword ptr [rsi-0x18]
| 2 | | | 0.2 | 0.8 | 1.0 | | | | | mov qword ptr [rbp+rdx*8-0x2040], r9 # wrong, HSW/SKL can keep indexed stores fused
| 1* | | | | | | | | | | mov edx, eax # gcc code-gen isn't great, but not as bad as in the SIMD loop. No extra movzx, but not taking advantage of AL/AH
| 1 | 0.4 | | | | | | 0.6 | | | shr eax, 0x18
| 1 | 0.8 | | | | | | 0.2 | | | shr edx, 0x10
| 1 | | 0.6 | | | | 0.3 | | | | add rax, 0x300
| 1* | | | | | | | | | | movzx edx, dl
| 1 | 0.2 | 0.1 | | | | 0.5 | 0.2 | | | add rdx, 0x200
| 1 | | | 0.5 0.5 | 0.5 0.5 | | | | | | mov r9, qword ptr [rbp+rdx*8-0x2040]
| 2^ | | 0.6 | 0.5 0.5 | 0.5 0.5 | | 0.3 | 0.1 | | | xor r9, qword ptr [rsi-0x10]
| 2 | | | 0.5 | 0.5 | 1.0 | | | | | mov qword ptr [rbp+rdx*8-0x2040], r9 # wrong, HSW/SKL can keep indexed stores fused
| 1 | | | 0.5 0.5 | 0.5 0.5 | | | | | | mov rdx, qword ptr [rbp+rax*8-0x2040]
| 2^ | | | 0.5 0.5 | 0.5 0.5 | | 0.6 | 0.4 | | | xor rdx, qword ptr [rsi-0x8]
| 2 | | | 0.5 | 0.5 | 1.0 | | | | | mov qword ptr [rbp+rax*8-0x2040], rdx # wrong, HSW/SKL can keep indexed stores fused
| 1 | 0.6 | | | | | | 0.4 | | | cmp r8, rdi
| 0F | | | | | | | | | | jnz 0xffffffffffffff75
Total Num Of Uops: 33
该循环比 IACA 的计数短 4 个融合域 uops，因为它不知道只有 SnB/IvB 非层压索引存储。HSW/SKL 没有。但是，这样的商店仍然不能使用端口 7，因此对于 4 个元素，这不会比 ~6.5 个周期更好。
（顺便说一句，通过对 Indices[i] 的简单处理，使用 movzx 分别加载每个元素，您将获得 4 个元素的 8 个周期，使端口 2 和 3 饱和。即使 gcc 不会生成用于解包结构的吞吐量优化代码， 4 字节的加载 + 解包应该是通过减轻一些加载端口压力而获得的净胜利。）
清理循环：
AVX2 在这里非常出色：我们循环遍历直方图的最低切片，并在其他切片中进行 XOR。这个循环是 8 个前端 uop，在 Skylake 上有 4 个负载，应该以每 2 个时钟 1 个迭代器运行：
.L7:
vmovdqa ymm2, YMMWORD PTR [rax+4096]
vpxor ymm0, ymm2, YMMWORD PTR [rax+6144]
vmovdqa ymm3, YMMWORD PTR [rax]
vpxor ymm1, ymm3, YMMWORD PTR [rax+2048]
vpxor ymm0, ymm0, ymm1
vmovdqa YMMWORD PTR [rax], ymm0
add rax, 32
cmp rax, rdx
jne .L7
我试图通过在一个链中执行 XOR 来进一步减少 uop 计数，但是 gcc 坚持执行两次vmovdqa加载并且必须在vpxor没有内存操作数的情况下执行一次。（OoO exec 将隐藏这个 VPXOR 的小链/树的延迟，所以没关系。）
我将如何使用 AVX-512 的散射？是否有一些分散指令可以异或而不是覆盖？
不，您将使用收集来获取旧值，然后使用 SIMD XOR，然后将更新的元素分散回它们来自的位置。
为避免冲突，您可能希望out[8][256]每个向量元素都可以使用不同的表。（否则，如果Indices[i+0]和Indices[i+4]相等，您就会遇到问题，因为分散存储只会存储具有该索引的最高向量元素。
分散/收集指令需要一个基址寄存器，但您可以_mm256_setr_epi64(0, 256, 256*2, ...);在进行vpmovzxbq零扩展加载后简单地添加。
笔记
我使用 IACA2.3 进行标量分析，因为当您在一个文件中有多个标记时，IACA3.0 似乎删除了-mark选择要分析的循环的选项。在这种情况下，IACA3.0 没有修复 IACA2.3 对 SKL 管道的任何错误。
推荐阅读
- flutter - 当我使用 list.filled() 并初始化值时，为什么我仍然提示语法错误？
- visual-studio - Visual Studio 2019 - 添加 Web 引用问题
- postgresql - 如何将时间四舍五入到 postgreSQL 中特定的最近给定间隔？
- sql - 如何在多个日期行上透视表
- ignite - 在 Apache Ignite 中，如何控制在哪个节点上创建缓存
- c# - 如何使用 OpenTK 渲染四边形？
- reactjs - 如何在反应路由中使用上下文提供者（反应 17）
- sql - 语法错误：应为“（”或“，”或关键字 SELECT，但脚本结束
- php - 如何在 Web 应用程序中实现版本控制？
- java - 从 maven 插件中提取依赖项的许可证