c++ - 我可以将AVX/SSE用作AoS布局而不是SoA吗？

Question

我想加速一个简单的积分器，该积分器通过其位置和速度描述一组无质量的粒子。我不是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进行上述计算，或者可能我错过了内在函数？

Answer 1

有关某些链接，请参见sse标签Wiki，尤其是SIMD at Insomniac Games (GDC 2015)。循环很多粒子与循环游戏世界中的许多对象是完全一样的问题，因此提到了这种循环以及尝试使用AoS的问题。

实际上，您不需要数组结构； xy[i]和vxvy[i]之间的距离是编译时常数并不重要。 (尽管它可能潜在地节省了寄存器，并为另一个指针增加了一些循环开销。但是，严重的是，大多数人不使用巨型结构，如果在编译时不知道其大小，则使用单独分配的数组。不过，它们可能会将所有指针保留在一个结构中。)

您(或编译器)可以为AoS方法改组并获得标量加速，但是无论如何，如果您遍历每个粒子，那么您就将自己射在了脚上。您的float2 xy对仅以64位块的形式出现，因此您无法有效地使用128位存储。使用两倍多的64位存储很糟糕:您将损失的SSE仅为一半，即AVX的75％。最重要的是，您要花费额外的指令进行改组或加载以及存储。

移动数据所花费的成本比实际的乘/加成多或多，特别是对于吞吐量而不是延迟。即使具有完美的SoA布局，也不会成为瓶颈的FMA单元，而是加载/存储吞吐量或总指令吞吐量(CPU前端)而不会显着循环展开。最重要的是，添加任何开销都意味着您只是前端的瓶颈(或洗牌吞吐量)。

当它们与vxvy交错时，无论是否显式存储到vxvy，都无法解决包含xy的缓存行的问题，因此对于此类问题，您总是需要两倍于SoA的存储带宽。

使用AVX处理不良的AoS布局时，值得手动洗牌，然后进行256b存储，该存储将存储新的xy值，并使用您先前加载的值重写vxvx，但是the compiler isn't allowed to do that when auto-vectorizing除非整个程序优化证明该代码是单线程。 C11和C++ 11内存模型一致认为，一个线程编写某些数组元素或struct成员，而其他线程编写其他元素则不是一场数据竞赛。当源仅读取那些元素时，不允许对vxvy成员进行非原子的读取-修改-写入。 (即，即使编译器重写了原始代码中的数据，也不允许编译器发明对不是由原始代码编写的内存位置/对象的写操作。)当然，如果您使用内在函数手动进行操作，则编译器可以做到这一点。如果需要，甚至particles[i].vxvy = particles[i].vxvy;也会授予编译器读取/随机播放/重写的许可。

实际上，编译器可以通过使用vmaskmovps进行掩码存储(而不是常规vmovups存储)来向量化此方式。它比常规商店要慢(Haswell / Skylake:需要p0，p1，p4和p23的4个融合域uops，而常规商店是需要p4和p237的单个微融合uop)。 Normally you want to avoid it，但是当编译器需要避免重写vxvy字节时，使用它进行矢量化可能仍然比单独的64位存储更好。特别是对于AVX512，带掩码的存储将允许使用512b(64字节) vector 自动矢量化，该 vector 一次存储4对xy对。 (对于SoA格式，则为8)。

我检查了gcc和ICC如何自动矢量化您的第一个版本，其中xy仍在AoS中，但布局与vxvy匹配，因此它可以使用纯垂直SIMD ops自动矢量化。 (source + asm output on the Godbolt compiler explorer)。 gcc正常，用一个vfmadd213ps指令进行循环。 ICC似乎对float2结构感到困惑，并且(我认为)实际上是在乘法/加法之前改组去交错，然后在之后再交错! (我不让ICC使用AVX或AVX512，因为更长的 vector 意味着更多的改组，因此甚至很难看清它在做什么。)这是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)。