我可以使用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,尤其是Insomniac Games(GDC 2015)的SIMD.循环很多粒子与在游戏世界中循环很多对象的问题完全相同,所以提到了这种循环,以及尝试使用AoS的问题.

你根本不需要数组结构; 并且xy[i]和之间的距离vxvy[i]是编译时常量并不重要.(但它可能会保存一个寄存器,并为另一个指针增量提供一些循环开销.但严重的是,大多数人不使用巨型结构,如果在编译时不知道大小,它们会使用单独分配的数组.但是,它们可能会将所有指针保留在结构中.)

对于你的AoS方法,你(或编译器)可以随机播放并获得超过标量的加速,但是如果你无论如何都要在每个粒子上循环,那么你就是通过这样做来拍摄自己.您的float2 xy对只有64位块,因此您无法使用128位存储.使用两倍多的64位存储很糟糕:你失去了SSE功耗的一半,或者是AVX功率的75%.最重要的是,你需要花费额外的指令进行改组或加载以及存储.

移动数据的成本与实际乘法/加法相同或更多,特别是对于吞吐量而非延迟.即使具有完美的SoA布局,也不是​​FMA单元是瓶颈,它将是加载/存储吞吐量,或者没有显着的循环展开的总指令吞吐量(CPU前端).在此之上添加任何开销意味着您只是前端的瓶颈(或随机吞吐量).

没有办法弄乱缓存行包含vxvy是否存储到vxvy显式与否,当它们交错时xy 因此,对于像这样的问题,你总是需要两倍的存储带宽与SoA相比.

使用AVX来处理糟糕的AoS布局,手动随机播放然后执行存储新xy值的256b存储并vxvx使用之前加载的值重写它是值得的,但是在自动向量化时不允许编译器执行此操作除非整个程序优化证明此代码是单线程的.C11和C++ 11内存模型同意,当一个线程编写某些数组元素或结构成员而其他线程正在编写其他元素时,它不是数据竞争.vxvy当源只读取这些元素时,不允许对成员进行非原子读 - 修改 - 写.(即,不允许编译器发明写入未由原始代码写入的内存位置/对象,即使它正在重写原来的数据.)当然,如果您使用内在函数手动执行,编译器可以执行此操作.particles[i].vxvy = particles[i].vxvy;如果需要,甚至可以给编译器许可证进行读/重写/重写.

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

我查看了gcc和ICC如何自动矢量化你的第一个版本,它xy仍然在AoS中,但在匹配的布局中vxvy,所以它可以使用纯垂直SIMD操作自动矢量化.(Godbolt编译器资源管理器上的source + asm输出).gcc没问题,用一条vfmadd213ps指令做一个循环.ICC似乎对float2结构感到困惑,并且(我认为)实际上是在乘法/加法之前进行去交错,然后重新交错!(我没有让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个融合域uop).