我可以将 AVX/SSE 与 swizzling 一起用于 AoS 布局而不是 SoA 吗？
Can I use AVX/SSE with swizzling for an AoS layout, instead of SoA?
我想加速一个简单的积分器，它通过位置和速度描述一组无质量粒子。我不是 SSE/AVX 专家，但我发现 SIMD 扩展在这里可以产生什么很有趣。
许多论文建议使用数组结构：
struct {
static float2 xy[OUGHTA_BE_ENOUGH];
static float2 vxvy[OUGHTA_BE_ENOUGH];
} Particles;
// in main loop:
Particles.xy[i] += time_delta * Particles.vxvy[i];
但是，对于许多应用程序而言，相反的方法将是有益的：
struct {
float2 xy;
float2 vxvy;
} Particle;
// in main loop:
particles[i].xy += time_delta * particles[i].vxvy;
虽然我模糊地理解要搜索什么来矢量化数组结构版本，但我怀疑是否有任何方法可以将 SIMD 与结构数组版本一起使用，因为字段访问或“调配”。
是否有任何技术可以使用 SIMD 进行上述计算，或者可能有我遗漏的内在函数？
论文中建议的方法可能更好，因为假设您使用 SSE 进行编译，编译器可以同时处理 4 个浮点运算。因此，循环可能会扩展为同时处理两个粒子。使用 gcc 和 O3 它确实发生了，只要看看反汇编就可以了。
使用您建议的方法，这仍然是可能的，但需要大量开销。
如果您担心保持单一线性内存访问，以我对现代处理器的理解，这不会有什么不同。处理器可以处理与缓存相关的多个线性内存访问。
现在，由于您正在模拟一组粒子，我建议对多个粒子进行矢量化。
struct Pack
{
floatN x; // N = 2, 4, 8 depending on your vectorization target
floatN y;
};
这样你就可以同时处理多个粒子，并且更容易编写矢量化。
f(Pack& pos, Pack speed) {
pos.x = add(pos.x, speed.x);
pos.y = add(pos.y, speed.y);
}
缺点是如果计算同时涉及两个粒子，它可能不适用：particleASpeed += force(particleBPos, particleAPos)。
此外，考虑使用 OpenCL 或 Cuda 进行此类计算，并考虑 运行 在 GPU 上使用它们。
最后，尽管它可能来得更早，但请记住衡量您在前后尝试优化的内容。我的建议充其量只是有根据的猜测，可能不会根据您的实际实施和问题的大小而改变任何内容。
见sse tag wiki for some links, especially SIMD at Insomniac Games (GDC 2015)。在很多粒子上循环与在游戏世界中在很多对象上循环是完全相同的问题，所以提到了这种循环，以及尝试使用 AoS 的问题。
您实际上不需要数组结构； xy[i] 和 vxvy[i] 之间的距离是编译时常量并不重要。 （不过，它可能会节省一个寄存器，并为另一个指针增量节省一些循环开销。但说真的，大多数人不使用巨型结构，如果在编译时不知道大小，他们会使用单独分配的数组。不过，他们可能会将所有 指针 保留在一个结构中。）
你（或编译器）可以为你的 AoS 方法洗牌并获得超过标量的加速，但如果你无论如何都要遍历每个粒子，你就是在开枪通过这样做在脚上。您的 float2 xy 对仅以 64 位块的形式出现，因此您不能有效地使用 128 位存储。使用两倍多的 64 位存储很糟糕：您失去了 SSE 的一半功能，或者 AVX 的 75% 的功能。除此之外，您还要花费额外的指令来改组或加载以及存储。
移动数据的成本与实际 multiply/add 一样多或更多，尤其是吞吐量而不是延迟。即使有一个完美的 SoA 布局，瓶颈也不会是 FMA 单元，而是 load/store 吞吐量，或没有显着循环展开的总指令吞吐量（CPU 前端）。在此之上添加任何开销意味着您只是在前端（或洗牌吞吐量）出现瓶颈。
无论您是否显式存储到 vxvy，当它们与 xy[= 交错时，都无法避免弄脏包含 vxvy 的缓存行68=] 因此，对于这样的问题，您总是需要两倍于 SoA 的存储带宽。
使用 AVX 解决你糟糕的 AoS 布局，手动洗牌然后做一个 256b 存储来存储新的 xy 值并用你之前加载的值重写 vxvx 是值得的, 但 除非全程序优化证明这段代码是单线程的。 C11 和 C++11 内存模型一致认为，当其他线程正在写入 other 元素时，一个线程写入一些数组元素或结构成员不是数据竞争。当源仅读取这些元素时，不允许对 vxvy 成员进行非原子读-修改-写。 （即不允许编译器发明写入不是由原始代码写入的内存位置/对象，即使它正在重写最初存在的数据。）当然，如果您使用内在函数手动执行，编译器可以做到.如果需要，甚至 particles[i].vxvy = particles[i].vxvy; 可能会授予编译器读取/随机播放/重写的许可。
实际上，编译器可以通过使用 vmaskmovps 进行掩码存储而不是常规 vmovups 存储 以这种方式矢量化。它比常规存储慢（Haswell/Skylake：需要 p0、p1、p4 和 p23 的 4 个融合域微指令，而常规存储是需要 p4 和 p237 的单个微融合微指令）。 ，但当编译器 需要 避免重写 vxvy 字节时，用它自动矢量化它可能仍然比单独的 64 位存储更好。特别是对于 AVX512，屏蔽存储将允许使用 512b（64 字节）向量进行自动向量化，一次存储 4 xy 对。 （而不是 8 个 SoA 格式）。
我检查了 gcc 和 ICC 如何自动矢量化您的第一个版本，其中 xy 仍在 AoS 中，但布局与 vxvy 匹配，因此它可以自动矢量化纯垂直 SIMD 操作。 (source + asm output on the Godbolt compiler explorer)。 gcc 没问题，用一条 vfmadd213ps 指令进行循环。 ICC 似乎对 float2 结构感到困惑，并且（我认为）实际上在 multiply/add 之前洗牌以去交错，然后在之后重新交错！ （我没有让 ICC 使用 AVX 或 AVX512，因为更长的矢量意味着更多的混洗，所以更难看出它在做什么。）这是 ICC 自动矢量化比 gcc 差的罕见情况之一。
gcc 和 ICC 都无法自动向量化 update_aos。这是我手动对其进行矢量化的方式（对于 AVX + FMA）：
// struct definitions and float2 operator*(float scalar, const float2 &f2)
// included in the Godbolt link, see above.
#include <immintrin.h>
void update_aos_manual(Particle *particles, size_t size, float time_delta_scalar)
{
__m256 time_delta = _mm256_set1_ps(time_delta_scalar);
// note: compiler can't prove this loop isn't infinite. (because i+=2 could wrap without making the condition false)
for(size_t i=0 ; i<size ; i+=2) {
float *ptr = (float*)(particles + i);
__m256 p = _mm256_load_ps(ptr); // xy0 vxvx0 | xy1 vxvy1
__m256 vx0 = _mm256_unpackhi_ps(p, _mm256_setzero_ps()); // vx0 0 | vx1 0
p = _mm256_fmadd_ps(time_delta, vx0, p); // p = td*vx0 + p
_mm256_store_ps(ptr, p);
//particles[i].xy += time_delta * particles[i].vxvy;
//particles[i].vxvy += 0.0f * particles[i].vxvy;
}
}
使用 gcc 和 ICC，这会编译成一个内部循环，如
## gcc7.2 -O3 -march=haswell
# various scalar setup for the loop:
vbroadcastss ymm0, xmm0 # ymm0 set1(time_delta_scalar)
vxorps xmm3, xmm3, xmm3 # ymm3 = setzero_ps
.L27:
vmovaps ymm2, YMMWORD PTR [rdi] # load 2 struct Particle
add rdi, 32 # ptr+=2 (16 bytes per element)
vunpckhps ymm1, ymm2, ymm3 # unpack high half of each lane with zero
vfmadd132ps ymm1, ymm2, ymm0 # xy += delta*vxvy; vxvy += delta*0
vmovaps YMMWORD PTR [rdi-32], ymm1 # store back into the array
cmp rdi, rax
jne .L27
这浪费了一半的存储带宽（不可避免）和一半的 FMA 吞吐量，但我认为您无法做得更好。好吧，展开会有所帮助，但是改组/取消改组和使用更少的 FMA 可能仍然是前端的瓶颈。通过展开，您可以在 Haswell/Skylake.
上以每个时钟几乎一个 32B 存储（每个时钟 4 个融合域 uops）将其存储到 运行
我想加速一个简单的积分器，它通过位置和速度描述一组无质量粒子。我不是 SSE/AVX 专家，但我发现 SIMD 扩展在这里可以产生什么很有趣。
许多论文建议使用数组结构：
struct {
static float2 xy[OUGHTA_BE_ENOUGH];
static float2 vxvy[OUGHTA_BE_ENOUGH];
} Particles;
// in main loop:
Particles.xy[i] += time_delta * Particles.vxvy[i];
但是，对于许多应用程序而言，相反的方法将是有益的：
struct {
float2 xy;
float2 vxvy;
} Particle;
// in main loop:
particles[i].xy += time_delta * particles[i].vxvy;
虽然我模糊地理解要搜索什么来矢量化数组结构版本，但我怀疑是否有任何方法可以将 SIMD 与结构数组版本一起使用，因为字段访问或“调配”。
是否有任何技术可以使用 SIMD 进行上述计算，或者可能有我遗漏的内在函数？
论文中建议的方法可能更好，因为假设您使用 SSE 进行编译，编译器可以同时处理 4 个浮点运算。因此，循环可能会扩展为同时处理两个粒子。使用 gcc 和 O3 它确实发生了，只要看看反汇编就可以了。
使用您建议的方法，这仍然是可能的，但需要大量开销。
如果您担心保持单一线性内存访问，以我对现代处理器的理解，这不会有什么不同。处理器可以处理与缓存相关的多个线性内存访问。
现在，由于您正在模拟一组粒子，我建议对多个粒子进行矢量化。
struct Pack
{
floatN x; // N = 2, 4, 8 depending on your vectorization target
floatN y;
};
这样你就可以同时处理多个粒子，并且更容易编写矢量化。
f(Pack& pos, Pack speed) {
pos.x = add(pos.x, speed.x);
pos.y = add(pos.y, speed.y);
}
缺点是如果计算同时涉及两个粒子，它可能不适用：particleASpeed += force(particleBPos, particleAPos)。
此外，考虑使用 OpenCL 或 Cuda 进行此类计算，并考虑 运行 在 GPU 上使用它们。
最后，尽管它可能来得更早，但请记住衡量您在前后尝试优化的内容。我的建议充其量只是有根据的猜测，可能不会根据您的实际实施和问题的大小而改变任何内容。
见sse tag wiki for some links, especially SIMD at Insomniac Games (GDC 2015)。在很多粒子上循环与在游戏世界中在很多对象上循环是完全相同的问题，所以提到了这种循环，以及尝试使用 AoS 的问题。
您实际上不需要数组结构； xy[i] 和 vxvy[i] 之间的距离是编译时常量并不重要。 （不过，它可能会节省一个寄存器，并为另一个指针增量节省一些循环开销。但说真的，大多数人不使用巨型结构，如果在编译时不知道大小，他们会使用单独分配的数组。不过，他们可能会将所有 指针 保留在一个结构中。）
你（或编译器）可以为你的 AoS 方法洗牌并获得超过标量的加速，但如果你无论如何都要遍历每个粒子，你就是在开枪通过这样做在脚上。您的 float2 xy 对仅以 64 位块的形式出现，因此您不能有效地使用 128 位存储。使用两倍多的 64 位存储很糟糕：您失去了 SSE 的一半功能，或者 AVX 的 75% 的功能。除此之外，您还要花费额外的指令来改组或加载以及存储。
移动数据的成本与实际 multiply/add 一样多或更多，尤其是吞吐量而不是延迟。即使有一个完美的 SoA 布局，瓶颈也不会是 FMA 单元，而是 load/store 吞吐量，或没有显着循环展开的总指令吞吐量（CPU 前端）。在此之上添加任何开销意味着您只是在前端（或洗牌吞吐量）出现瓶颈。
无论您是否显式存储到 vxvy，当它们与 xy[= 交错时，都无法避免弄脏包含 vxvy 的缓存行68=] 因此，对于这样的问题，您总是需要两倍于 SoA 的存储带宽。
使用 AVX 解决你糟糕的 AoS 布局，手动洗牌然后做一个 256b 存储来存储新的 xy 值并用你之前加载的值重写 vxvx 是值得的, 但 vxvy 成员进行非原子读-修改-写。 （即不允许编译器发明写入不是由原始代码写入的内存位置/对象，即使它正在重写最初存在的数据。）当然，如果您使用内在函数手动执行，编译器可以做到.如果需要，甚至 particles[i].vxvy = particles[i].vxvy; 可能会授予编译器读取/随机播放/重写的许可。
实际上，编译器可以通过使用 vmaskmovps 进行掩码存储而不是常规 vmovups 存储 以这种方式矢量化。它比常规存储慢（Haswell/Skylake：需要 p0、p1、p4 和 p23 的 4 个融合域微指令，而常规存储是需要 p4 和 p237 的单个微融合微指令）。 vxvy 字节时，用它自动矢量化它可能仍然比单独的 64 位存储更好。特别是对于 AVX512，屏蔽存储将允许使用 512b（64 字节）向量进行自动向量化，一次存储 4 xy 对。 （而不是 8 个 SoA 格式）。
我检查了 gcc 和 ICC 如何自动矢量化您的第一个版本，其中 xy 仍在 AoS 中，但布局与 vxvy 匹配，因此它可以自动矢量化纯垂直 SIMD 操作。 (source + asm output on the Godbolt compiler explorer)。 gcc 没问题，用一条 vfmadd213ps 指令进行循环。 ICC 似乎对 float2 结构感到困惑，并且（我认为）实际上在 multiply/add 之前洗牌以去交错，然后在之后重新交错！ （我没有让 ICC 使用 AVX 或 AVX512，因为更长的矢量意味着更多的混洗，所以更难看出它在做什么。）这是 ICC 自动矢量化比 gcc 差的罕见情况之一。
gcc 和 ICC 都无法自动向量化 update_aos。这是我手动对其进行矢量化的方式（对于 AVX + FMA）：
// struct definitions and float2 operator*(float scalar, const float2 &f2)
// included in the Godbolt link, see above.
#include <immintrin.h>
void update_aos_manual(Particle *particles, size_t size, float time_delta_scalar)
{
__m256 time_delta = _mm256_set1_ps(time_delta_scalar);
// note: compiler can't prove this loop isn't infinite. (because i+=2 could wrap without making the condition false)
for(size_t i=0 ; i<size ; i+=2) {
float *ptr = (float*)(particles + i);
__m256 p = _mm256_load_ps(ptr); // xy0 vxvx0 | xy1 vxvy1
__m256 vx0 = _mm256_unpackhi_ps(p, _mm256_setzero_ps()); // vx0 0 | vx1 0
p = _mm256_fmadd_ps(time_delta, vx0, p); // p = td*vx0 + p
_mm256_store_ps(ptr, p);
//particles[i].xy += time_delta * particles[i].vxvy;
//particles[i].vxvy += 0.0f * particles[i].vxvy;
}
}
使用 gcc 和 ICC，这会编译成一个内部循环，如
## gcc7.2 -O3 -march=haswell
# various scalar setup for the loop:
vbroadcastss ymm0, xmm0 # ymm0 set1(time_delta_scalar)
vxorps xmm3, xmm3, xmm3 # ymm3 = setzero_ps
.L27:
vmovaps ymm2, YMMWORD PTR [rdi] # load 2 struct Particle
add rdi, 32 # ptr+=2 (16 bytes per element)
vunpckhps ymm1, ymm2, ymm3 # unpack high half of each lane with zero
vfmadd132ps ymm1, ymm2, ymm0 # xy += delta*vxvy; vxvy += delta*0
vmovaps YMMWORD PTR [rdi-32], ymm1 # store back into the array
cmp rdi, rax
jne .L27
这浪费了一半的存储带宽（不可避免）和一半的 FMA 吞吐量，但我认为您无法做得更好。好吧，展开会有所帮助，但是改组/取消改组和使用更少的 FMA 可能仍然是前端的瓶颈。通过展开，您可以在 Haswell/Skylake.
上以每个时钟几乎一个 32B 存储（每个时钟 4 个融合域 uops）将其存储到 运行