使用 AVX2 指令选择性地异或列表元素
selectively xor-ing elements of a list with AVX2 instructions
我想用 AVX2 指令加速以下操作，但我找不到方法。
我得到了一个 uint64_t 的大数组 uint64_t data[100000] 和一个 unsigned char indices[100000] 字节的数组。我想输出一个数组 uint64_t Out[256]，其中第 i 个值是所有 data[j] 的异或，使得 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];
}
您可以根据索引[i]对数据进行排序...这应该需要 O(N*log2(N))，但可以并行化。
然后对排序后的数据进行累积异或——也可以并行化。
那么就是计算Out[i] = CumXor(j) ^ Out[i-1];
的事情了
基于 Haswell/Skylake 的静态分析，我想出了一个版本，当由海合会。大尺寸的平均值，不包括组合 Out[] 的多个副本的时间，并假设索引的随机分布不会导致任何 store/reload 深度链 运行ning 足够长重要。
如果您关心 Ryzen 或 Excavator（其他 2 个主流 AVX2 微架构），请 IDK。
我没有手工仔细分析过，但是 IACA 对于 HSW/SKL 是错误的，并且认为有些指令实际上没有微熔断（在 i7 上测试过） -6700k with perf counters)，所以它认为前端瓶颈比实际情况更严重。例如movhps load+merge micro-fuses，但 IACA 认为它甚至没有简单的寻址模式。
我们应该忽略任何缓存未命中，因为 uint64_t Out[4][256] 只有 8kiB。因此，我们的缓存占用空间仅为最新 CPU 上 L1d 大小的 1/4，即使超线程在两个逻辑线程之间共享 L1d，也应该基本没问题。遍历 data[] 和 Indices[] 应该可以很好地预取，并且希望不会驱逐 Out[] 太多。因此，静态分析很有可能在某种程度上是准确的，而且它比仔细的微基准测试更快，更重要的是，它可以准确地告诉您瓶颈是什么。
当然，我们严重依赖乱序执行，不完善的调度或其他意想不到的瓶颈很容易发生。不过，如果我没有得到报酬，我真的不想进行微基准测试。
Can we implement this more efficiently with AVX2 instructions?
这基本上是一个直方图问题。使用多个 table 并在最后组合的常用直方图优化适用 。 SIMD XOR 对于最后合并很有用（只要您使用 Out[4][256]，而不是 Out[256][4]。后者还需要通过 8*4 而不是缩放来使索引变慢8（可以在缩放索引寻址模式下使用单个 LEA 完成））。
但与普通直方图不同的是，您要对内存中的一些数据进行异或运算，而不是将常数 1 添加。因此，代码必须将 data[i] 加载到注册为 xor 的来源。 （或者加载，然后 xor reg, data[i]/store）。这比直方图的总内存操作还要多。
我们从"manual"gather/scatter提前进入SIMD向量（使用movq/movhpsloads/stores），允许我们将 SIMD 用于 data[i] 加载和 XOR。这减少了加载操作的总数，从而降低了加载端口压力，而无需花费额外的前端带宽。
手动收集到 256 位向量可能不值得额外的洗牌（一个额外的 vinserti128 / vextracti128 只是为了我们可以将 2 个内存源 vpxor 组合成一个 256 位的）。 128位向量应该不错。前端吞吐量也是一个主要问题，因为（在 Intel SnB 系列 CPU 上）您希望避免存储的索引寻址模式。 gcc 使用 lea 指令计算寄存器中的地址，而不是使用索引 loads/stores。 -march=skylake 的 clang / LLVM 决定不这样做，在这种情况下这是一个错误的决定，因为端口 2 / 端口 3 上的循环瓶颈，并且花费额外的 ALU 微指令来允许存储地址微指令使用端口 7 是一个胜利.但是，如果您 没有 在 p23 上遇到瓶颈，那么花费额外的 uops 来避免索引存储是不好的。 （还有 in cases where the can stay micro-fused，绝对不仅仅是为了避免索引加载；愚蠢的 gcc）。也许 gcc 和 LLVM 的寻址模式成本模型不是很准确，或者它们没有足够详细地模拟管道以找出前端与特定端口何时出现循环瓶颈。
寻址模式的选择和其他 asm 代码生成选择对于在 SnB 系列上实现最佳性能至关重要。但是用 C 语言编写你无法控制它；您主要受编译器的支配，除非您可以调整源代码以使其做出不同的选择。例如gcc vs. clang 在这里有很大的不同。
在 SnB 系列上，movhps 负载需要端口 5 用于 shuffle/blend（尽管它确实微融合到一个 uop），但是 movhps 存储是一个纯没有 ALU 微指令的存储。所以它是收支平衡的，让我们对两个数据元素使用一个 SIMD 加载/XOR。
使用 AVX，ALU 微指令允许未对齐的内存源操作数，因此我们不需要为 data[] 要求对齐。但是 Intel HSW/SKL 可以保持索引寻址模式与 pxor 微融合，而不是 vpxor。因此，在不启用 AVX 的情况下进行编译可能 更好 ，允许编译器使用索引寻址模式而不是递增单独的指针。 （或者如果编译器不知道这一点并使用索引寻址模式，则使其更快。）TL:DR：可能最好要求 16 字节对齐 data[] 并编译该函数禁用 AVX， 以获得更好的宏融合。 （但随后我们错过了 256 位 SIMD，用于在最后组合 Out 切片，除非我们将其放入使用 AVX 或 AVX2 编译的不同函数中）
避免未对齐的加载也将避免任何缓存行拆分，这不会花费额外的 uops，但我们可能接近 L1d 吞吐量限制的瓶颈，而不仅仅是 load/store 执行单元吞吐量限制。
我还查看了一次加载 4 个索引并使用 ALU 指令解包。例如用 memcpy 变成 struct { uint8_t idx[4]; } idx;。但是 gcc 会生成多个用于解包的浪费指令。太糟糕了 x86 没有像 ARM ubfx 或 尤其是 PowerPC rlwinm 这样好的位域指令（这可能会使结果自由左移，所以如果 x86如果那样的话，静态 Out 可以在非 PIC 代码中使用 base+disp32 寻址模式。）
如果我们使用标量异或，使用 shift / movzx 从 AL/AH 解压双字是成功的，但当我们对 data[] 使用 SIMD 并花费前端时，它看起来不是lea 指令上的 -end 吞吐量，以允许存储地址微指令到端口 7 上的 运行。这使我们成为前端瓶颈而不是 port2/3 瓶颈，因此使用 4x movzx根据静态分析，从内存加载看起来最好。如果您花时间手动编辑 asm，那么这两种方式都值得进行基准测试。 （gcc 生成的带有额外 uops 的 asm 很糟糕，包括在右移 24 位后完全冗余的 movzx，高位已经为零。）
密码
（在 the Godbolt compiler explorer 上查看它以及标量版本）：
#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，对索引和 data[] 使用相同的计数器，这样可以节省 add.
clang 不使用 LEA，瓶颈在每 7 个周期 4 is，因为 none 的存储 uops 可以 运行 在端口 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 个融合域微指令，因为它不知道只有 SnB/IvB 未层压索引存储。 HSW/SKL 不要。但是，这样的商店仍然不能使用端口 7，因此对于 4 个元素，这不会比 ~6.5 个周期更好。
（顺便说一句，通过简单处理索引 [i]，用 movzx 分别加载每个索引，你会得到 4 个元素的 8 个周期，使端口 2 和 3 饱和。即使 gcc 不生成吞吐量最佳代码对于解包结构，4 字节加载 + 解包应该是通过减轻一些加载端口压力而获得的净胜利。）
清理循环:
AVX2 在这里真的很出色：我们遍历直方图的最低切片，并在其他切片中进行 XOR。这个循环是 8 个前端微指令，在 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
我试图通过在一个链中进行异或来进一步减少 uop 计数，但是 gcc 坚持要进行两次 vmovdqa 加载并且必须在没有内存操作数的情况下进行一次 vpxor。 （OoO exec 将隐藏 VPXOR 的这个小链/树的延迟，所以没关系。）
How would I use the scattering with AVX-512? Is there some scatters instruction that xors instead of overwrites?
不，您会使用收集来获取旧值，然后进行 SIMD 异或，然后将更新的元素分散回它们来自的位置。
为避免冲突，您可能需要 out[8][256]，这样每个向量元素都可以使用不同的 table。 （否则，如果 Indices[i+0] 和 Indices[i+4] 相等，就会出现问题，因为散点存储只会存储具有该索引的最高向量元素。
Scatter/gather 指令需要一个基址寄存器，但您可以在执行 vpmovzxbq 零扩展加载后简单地添加 _mm256_setr_epi64(0, 256, 256*2, ...);。
备注
我使用 IACA2.3 进行标量分析，因为 IACA3.0 似乎已经删除了 -mark 选项来选择当一个文件中有多个标记时要分析哪个循环。 IACA3.0 没有修复 IACA2.3 在这种情况下关于 SKL 管道的任何错误。
我想用 AVX2 指令加速以下操作，但我找不到方法。
我得到了一个 uint64_t 的大数组 uint64_t data[100000] 和一个 unsigned char indices[100000] 字节的数组。我想输出一个数组 uint64_t Out[256]，其中第 i 个值是所有 data[j] 的异或，使得 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];
}
您可以根据索引[i]对数据进行排序...这应该需要 O(N*log2(N))，但可以并行化。
然后对排序后的数据进行累积异或——也可以并行化。
那么就是计算Out[i] = CumXor(j) ^ Out[i-1];
基于 Haswell/Skylake 的静态分析，我想出了一个版本，当由海合会。大尺寸的平均值，不包括组合 Out[] 的多个副本的时间，并假设索引的随机分布不会导致任何 store/reload 深度链 运行ning 足够长重要。
如果您关心 Ryzen 或 Excavator（其他 2 个主流 AVX2 微架构），请 IDK。
我没有手工仔细分析过，但是 IACA 对于 HSW/SKL 是错误的，并且认为有些指令实际上没有微熔断（在 i7 上测试过） -6700k with perf counters)，所以它认为前端瓶颈比实际情况更严重。例如movhps load+merge micro-fuses，但 IACA 认为它甚至没有简单的寻址模式。
我们应该忽略任何缓存未命中，因为 uint64_t Out[4][256] 只有 8kiB。因此，我们的缓存占用空间仅为最新 CPU 上 L1d 大小的 1/4，即使超线程在两个逻辑线程之间共享 L1d，也应该基本没问题。遍历 data[] 和 Indices[] 应该可以很好地预取，并且希望不会驱逐 Out[] 太多。因此，静态分析很有可能在某种程度上是准确的，而且它比仔细的微基准测试更快，更重要的是，它可以准确地告诉您瓶颈是什么。
当然，我们严重依赖乱序执行，不完善的调度或其他意想不到的瓶颈很容易发生。不过，如果我没有得到报酬，我真的不想进行微基准测试。
Can we implement this more efficiently with AVX2 instructions?
这基本上是一个直方图问题。使用多个 table 并在最后组合的常用直方图优化适用 。 SIMD XOR 对于最后合并很有用（只要您使用 Out[4][256]，而不是 Out[256][4]。后者还需要通过 8*4 而不是缩放来使索引变慢8（可以在缩放索引寻址模式下使用单个 LEA 完成））。
但与普通直方图不同的是，您要对内存中的一些数据进行异或运算，而不是将常数 1 添加。因此，代码必须将 data[i] 加载到注册为 xor 的来源。 （或者加载，然后 xor reg, data[i]/store）。这比直方图的总内存操作还要多。
我们从"manual"gather/scatter提前进入SIMD向量（使用movq/movhpsloads/stores），允许我们将 SIMD 用于 data[i] 加载和 XOR。这减少了加载操作的总数，从而降低了加载端口压力，而无需花费额外的前端带宽。
手动收集到 256 位向量可能不值得额外的洗牌（一个额外的 vinserti128 / vextracti128 只是为了我们可以将 2 个内存源 vpxor 组合成一个 256 位的）。 128位向量应该不错。前端吞吐量也是一个主要问题，因为（在 Intel SnB 系列 CPU 上）您希望避免存储的索引寻址模式。 gcc 使用 lea 指令计算寄存器中的地址，而不是使用索引 loads/stores。 -march=skylake 的 clang / LLVM 决定不这样做，在这种情况下这是一个错误的决定，因为端口 2 / 端口 3 上的循环瓶颈，并且花费额外的 ALU 微指令来允许存储地址微指令使用端口 7 是一个胜利.但是，如果您 没有 在 p23 上遇到瓶颈，那么花费额外的 uops 来避免索引存储是不好的。 （还有 in cases where the can stay micro-fused，绝对不仅仅是为了避免索引加载；愚蠢的 gcc）。也许 gcc 和 LLVM 的寻址模式成本模型不是很准确，或者它们没有足够详细地模拟管道以找出前端与特定端口何时出现循环瓶颈。
寻址模式的选择和其他 asm 代码生成选择对于在 SnB 系列上实现最佳性能至关重要。但是用 C 语言编写你无法控制它；您主要受编译器的支配，除非您可以调整源代码以使其做出不同的选择。例如gcc vs. clang 在这里有很大的不同。
在 SnB 系列上，movhps 负载需要端口 5 用于 shuffle/blend（尽管它确实微融合到一个 uop），但是 movhps 存储是一个纯没有 ALU 微指令的存储。所以它是收支平衡的，让我们对两个数据元素使用一个 SIMD 加载/XOR。
使用 AVX，ALU 微指令允许未对齐的内存源操作数，因此我们不需要为 data[] 要求对齐。但是 Intel HSW/SKL 可以保持索引寻址模式与 pxor 微融合，而不是 vpxor。因此，在不启用 AVX 的情况下进行编译可能 更好 ，允许编译器使用索引寻址模式而不是递增单独的指针。 （或者如果编译器不知道这一点并使用索引寻址模式，则使其更快。）TL:DR：可能最好要求 16 字节对齐 data[] 并编译该函数禁用 AVX， 以获得更好的宏融合。 （但随后我们错过了 256 位 SIMD，用于在最后组合 Out 切片，除非我们将其放入使用 AVX 或 AVX2 编译的不同函数中）
避免未对齐的加载也将避免任何缓存行拆分，这不会花费额外的 uops，但我们可能接近 L1d 吞吐量限制的瓶颈，而不仅仅是 load/store 执行单元吞吐量限制。
我还查看了一次加载 4 个索引并使用 ALU 指令解包。例如用 memcpy 变成 struct { uint8_t idx[4]; } idx;。但是 gcc 会生成多个用于解包的浪费指令。太糟糕了 x86 没有像 ARM ubfx 或 尤其是 PowerPC rlwinm 这样好的位域指令（这可能会使结果自由左移，所以如果 x86如果那样的话，静态 Out 可以在非 PIC 代码中使用 base+disp32 寻址模式。）
如果我们使用标量异或，使用 shift / movzx 从 AL/AH 解压双字是成功的，但当我们对 data[] 使用 SIMD 并花费前端时，它看起来不是lea 指令上的 -end 吞吐量，以允许存储地址微指令到端口 7 上的 运行。这使我们成为前端瓶颈而不是 port2/3 瓶颈，因此使用 4x movzx根据静态分析，从内存加载看起来最好。如果您花时间手动编辑 asm，那么这两种方式都值得进行基准测试。 （gcc 生成的带有额外 uops 的 asm 很糟糕，包括在右移 24 位后完全冗余的 movzx，高位已经为零。）
密码
（在 the Godbolt compiler explorer 上查看它以及标量版本）：
#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，对索引和 data[] 使用相同的计数器，这样可以节省 add.
clang 不使用 LEA，瓶颈在每 7 个周期 4 is，因为 none 的存储 uops 可以 运行 在端口 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 个融合域微指令，因为它不知道只有 SnB/IvB 未层压索引存储。 HSW/SKL 不要。但是，这样的商店仍然不能使用端口 7，因此对于 4 个元素，这不会比 ~6.5 个周期更好。
（顺便说一句，通过简单处理索引 [i]，用 movzx 分别加载每个索引，你会得到 4 个元素的 8 个周期，使端口 2 和 3 饱和。即使 gcc 不生成吞吐量最佳代码对于解包结构，4 字节加载 + 解包应该是通过减轻一些加载端口压力而获得的净胜利。）
清理循环:
AVX2 在这里真的很出色：我们遍历直方图的最低切片，并在其他切片中进行 XOR。这个循环是 8 个前端微指令，在 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
我试图通过在一个链中进行异或来进一步减少 uop 计数，但是 gcc 坚持要进行两次 vmovdqa 加载并且必须在没有内存操作数的情况下进行一次 vpxor。 （OoO exec 将隐藏 VPXOR 的这个小链/树的延迟，所以没关系。）
How would I use the scattering with AVX-512? Is there some scatters instruction that xors instead of overwrites?
不，您会使用收集来获取旧值，然后进行 SIMD 异或，然后将更新的元素分散回它们来自的位置。
为避免冲突，您可能需要 out[8][256]，这样每个向量元素都可以使用不同的 table。 （否则，如果 Indices[i+0] 和 Indices[i+4] 相等，就会出现问题，因为散点存储只会存储具有该索引的最高向量元素。
Scatter/gather 指令需要一个基址寄存器，但您可以在执行 vpmovzxbq 零扩展加载后简单地添加 _mm256_setr_epi64(0, 256, 256*2, ...);。
备注
我使用 IACA2.3 进行标量分析，因为 IACA3.0 似乎已经删除了 -mark 选项来选择当一个文件中有多个标记时要分析哪个循环。 IACA3.0 没有修复 IACA2.3 在这种情况下关于 SKL 管道的任何错误。