快速将 2 个 double 数组交织成具有 2 个 float 和 1 个 int(循环不变)成员的结构数组,并使用 SIMD double->float 转换?

Fast interleave 2 double arrays into an array of structs with 2 float and 1 int (loop invariant) member, with SIMD double->float conversion?

我有一段代码是 x86 处理器上 C++ 应用程序 运行 的瓶颈,我们从两个数组中获取双精度值,转换为浮点数并存储在结构数组中。这是一个瓶颈的原因是它被调用时有非常大的循环,或者被调用了数千次。

有没有更快的方法来使用 SIMD Intrinsics 执行此复制和强制转换操作?我看过 this answer on faster memcpy 但没有提到演员表。

简单的 C++ 循环案例如下所示

        int _iNum;
        const unsigned int _uiDefaultOffset; // a constant 

        double * pInputValues1; // array of double values, count = _iNum;
        double * pInputValues2; 

        MyStruct * pOutput;    // array of outputs defined as
        // struct MyStruct 
        // { 
        //    float O1;
        //    float O2;
        //    unsigned int Offset;
        // };

        for (int i = 0; i < _iNum; ++i)
        {
            _pPoints[i].O1 = static_cast<float>(pInputValues1[i]);
            _pPoints[i].O2 = static_cast<float>(pInputValues2[i]);
            _pPoints[i].Offset = _uiDefaultOffset;
        }

注意: 结构格式为 [Float,Float,Int](24 字节)但我们可以(如果它有助于提高性能)添加一个额外的4 字节填充使其成为 32 字节。

这是对 SSE4.1 的尝试,没有 AVX(这样做比较棘手,到目前为止我想出了更多的混洗),并使用 12 字节/点格式:(未测试)

void test3(MyStruct * _pPoints, double * pInputValues1, double * pInputValues2) {
        // struct MyStruct 
        // { 
        //    float O1;
        //    float O2;
        //    unsigned int Offset;
        // };
    __m128 offset = _mm_castsi128_ps(_mm_cvtsi32_si128(_uiDefaultOffset));
    int i;
    for (i = 0; i < _iNum - 2; i += 2)
    {
        // read inputs and convert to float
        __m128d inA = _mm_loadu_pd(&pInputValues1[i]);
        __m128d inB = _mm_loadu_pd(&pInputValues2[i]);
        __m128 inAf = _mm_cvtpd_ps(inA);    // 0 0 A1 A0
        __m128 inBf = _mm_cvtpd_ps(inB);    // 0 0 B1 B0
        // shuffle B0 from place 0 to place 1, merge with offset
        __m128 tempA = _mm_shuffle_ps(inBf, offset, _MM_SHUFFLE(1, 0, 0, 0)); // 0 OF B0 B0
        // shuffle A1 from place 1 to place 0, merge with offset
        __m128 tempB = _mm_shuffle_ps(inAf, offset, _MM_SHUFFLE(1, 0, 1, 1)); // 0 OF A1 A1
        // replace B0 at place 0 with A0
        __m128 outA = _mm_blend_ps(tempA, inAf, 1);  // 0 OF B0 A0
        // replace A1 at place 1 with B1
        __m128 outB = _mm_blend_ps(tempB, inBf, 2);  // 0 OF B1 A1
        // store results
        _mm_storeu_ps(&_pPoints[i].O1, outA);
        _mm_storeu_ps(&_pPoints[i + 1].O1, outB);
    }
    // remaining iteration if _iNum is not even
    for (; i < _iNum; i++)
    {
        _pPoints[i].O1 = static_cast<float>(pInputValues1[i]);
        _pPoints[i].O2 = static_cast<float>(pInputValues2[i]);
        _pPoints[i].Offset = _uiDefaultOffset;
    }
}

这使用了shufps的能力来从两个不同的来源中进行选择来合并动态数据和常量偏移量,同样的洗牌也会移动每个需要移动的组中的浮点数。然后使用混合将单个浮标替换为已经在正确位置的另一个浮标。这需要 2 次随机播放和 2 次混合,还有一种方法是 3 次随机播放和零混合,但是在当前的英特尔处理器上,所有随机播放都转到 p5,而混合可以转到不同的端口。转换也已经使用了 p5,所以它被淹没了,使用混合应该更好。每次迭代仍然是 4 p5 微操作,因此处理每个项目至少需要 2 个周期,这不是很好。

主循环跳过最后的项目,这样它就不会越界写入,它会 slightly-overlapping 16 字节存储,在结构末尾写入 4 个字节。该部分被下一个存储的真实结果覆盖,但在数组末尾这样做可能很危险。

这个问题和memcpy不是很像。这一切都是关于优化 loop-invariant 整数成员的混洗 and/or 标量存储的交错。这使 SIMD 变得困难。

需要使用这种存储格式,其中 intfloat 成员交错排列吗?交错浮动已经够糟糕了。我假设稍后的一些代码将修改不同结构中的 ints,否则为每个元素复制它是没有意义的。

您能否以 4 个元素为一组工作,例如 struct { float a[4], b[4]; int i[4]; };,这样您就可以将 4x 连续 double 加载并转换为 4x float 并进行 128 位 SIMD 存储?访问单个 output-array "struct".

的所有 3 个成员时,您仍然会有一些空间局部性

无论如何,假设您的输出格式必须完全交错,我们不需要将其填充到 16 字节。 x86 CPU 可以有效地处理重叠的 16 字节存储以处理 12 字节结构,如@harold 的回答显示。 Cache-line 拆分的成本可能低于存储填充所需的额外内存带宽。

或者另一种策略是为花车和 int 使用单独的商店,因此您不需要重叠。我们可能可以将其优化到每 2 个周期 1 个结构的每个时钟周期 1 个存储瓶颈。 (或者略低,因为 IIRC cache-split 存储需要重放存储 uop,至少在 Intel CPU 上是这样。)我们也可以展开 4*12 = 3*16 字节并通过使用 SIMD 存储来保存 2 个整数存储浮动数据。 48 字节 = xyIx|yIxy|IxyI 有四个 I 元素作为四个结构的一部分,但它们足够接近我们可以用两个 _mm_storeu_si128( set1(offset) ) 内在函数存储所有 4 个。然后存储与之重叠的 xy 对。 16 字节边界用 | 标记。如果缓存行拆分是一个问题,我们可以为最后一个 对齐的向量做 2x 标量和一个 SIMD(如果输出数组是 16 字节对齐的)。或者在 Intel Haswell 和更高版本的 CPU 上,32 字节对齐的存储可能很好。


如果我们不小心,我们很容易在 Intel CPU 上出现洗牌吞吐量瓶颈,尤其是 Sandybridge-family(SnB 通过 Skylake/Coffee Lake),其中 FP 洗牌只能 运行在端口 5 上。这就是为什么我正在考虑 而不是 将所有内容混在一起,每个结构 1 个存储。

SIMD double->float conversion 花费 2 微指令:随机播放 + FP-math,因为浮点数是宽度的一半,指令将浮点数打包到向量寄存器的底部。

AVX在这里很有用,可以将4doubles转换成4floats的SIMD向量。

除此之外,我同意@harold 的观点,即 128 位向量可能是一个不错的选择。甚至 AVX2 也没有很好的 2 输入 lane-crossing 洗牌,AVX1 也非常有限。所以我们可以使用 256-bit -> 128-bit double->float 转换来提供基于 __m128.

的交错策略

vmovhps [mem], xmm 不会在 Intel CPU 上花费 shuffle uop,只是一个纯粹的存储,所以将 2 个向量混合在一起并将 [ B1 A1 B0 A0 ] 放入一个向量中,我们就可以使用两个 64 位没有任何额外洗牌的低半部和高半部的存储。

OTOH,@harold 的版本可能仍然更好。每 2 个结构 4 次洗牌可能比每 2 个结构 4 次存储更好,因为存储有时需要重放以进行缓存行拆分,但洗牌不需要。但是使用重叠存储技巧,每 2 个结构 3.5 或 3 个存储看起来是可行的。


或者这是另一个想法,它使用了上面的一些方法,但进行了一些混合以节省商店

我基本上是在编辑@harold 的代码以实现我在上面的文本中写的想法时想到的。在这里使用混合是减少存储和洗牌端口压力的好方法。

上面的一些想法仍然值得探索,尤其是做一个广泛的 set1(offset) 存储,然后将它与 64 位 vmovlps 存储重叠。 (通过 3x2 = 6 或 3x4 = 12 输出结构展开后,使其成为我们一次转换的 4 个双精度数的倍数。)12 * 12 = 144 字节,它是 16 的倍数,但不是 32 或 64 的倍数,所以我们至少可以随时知道我们相对于 16 字节边界的位置,但除非我们展开更多,否则不知道缓存行。 (可能留下更多需要清理的工作,并且膨胀 code-size。)

#include <immintrin.h>
#include <stddef.h>
#include <stdint.h>

struct f2u { 
  float O1, O2;
  unsigned int Offset;
};

// names with a leading _ at file scope are reserved for the implementation.
// fixed that portability problem for you.
static const unsigned uiDefaultOffset = 123;


// only requires AVX1
// ideally pA and pB should be 32-byte aligned.
// probably also dst 16-byte aligned is good.
void cvt_interleave_avx(f2u *__restrict dst, double *__restrict pA, double *__restrict pB, ptrdiff_t len)
{
    __m128 voffset = _mm_castsi128_ps(_mm_set1_epi32(uiDefaultOffset));

    // 48 bytes per iteration: 3x16 = 4x12
    ptrdiff_t i;
    for (i = 0; i < len - 3; i += 4)
    {
        // read inputs and convert to float
        __m256d inA = _mm256_loadu_pd(&pA[i]);
        __m256d inB = _mm256_loadu_pd(&pB[i]);
        __m128 inAf = _mm256_cvtpd_ps(inA);    // A3 A2 A1 A0
        __m128 inBf = _mm256_cvtpd_ps(inB);    // B3 B2 B1 B0

        // interleave to get XY pairs
        __m128 lo = _mm_unpacklo_ps(inAf, inBf); // B1 A1 B0 A0
        __m128 hi = _mm_unpackhi_ps(inAf, inBf); // B3 A3 B2 A2

        // blend integer into place
        __m128 out0 = _mm_blend_ps(lo, voffset, 1<<2);  // x OF B0 A0
        __m128 out2 = _mm_blend_ps(hi, voffset, 1<<2);  // x OF B2 A2

        // TODO: _mm_alignr_epi8 to create OF OF B1 A1 spending 1 more shuffle to save a store.

        // store results
        _mm_storeu_ps(&dst[i + 0].O1, out0);  // 16 bytes with blended integer
        _mm_storeh_pi((__m64*)&dst[i + 1].O1, lo);    // 8 bytes from top half of reg, partial overlap
        dst[i + 1].Offset = uiDefaultOffset;

        _mm_storeu_ps(&dst[i + 2].O1, out2);  // 16 bytes with blended integer
        _mm_storeh_pi((__m64*)&dst[i + 3].O1, hi);    // 8 bytes from top half of reg, partial overlap
        dst[i + 3].Offset = uiDefaultOffset;
    }

    // scalar cleanup for  if _iNum is not even
    for (; i < len; i++)
    {
        dst[i].O1 = static_cast<float>(pA[i]);
        dst[i].O2 = static_cast<float>(pB[i]);
        dst[i].Offset = uiDefaultOffset;
    }
}

gcc9.1 -O3 -march=skylake on Godbolt 将主循环编译为 front-end 的 19 fused-domain 微指令。 (vcvtpd2ps 指令都不能 micro-fuse 因为 GCC 没有做任何聪明的事情,比如相对于 pA 寻址 pB 来避免其中之一的索引寻址模式。所以他们'每 3 微秒:加载 + 转换 + 随机播放)

但无论如何它确实在 back-end 中的存储上出现瓶颈,即使每次迭代需要完整的 5 个周期才能从 4-wide front-end.

发出

每次迭代有 6 个存储(每 4 个结构),这将使它成为每 6 个周期最多 1 次迭代的瓶颈,瓶颈在 store-data port/execution 单元上。 (直到每个时钟可以做 2 次存储的 Ice Lake。)所以这实现了理论最佳情况下每 1.5 个周期 1 个结构, 与我之前为 overlapping-store 想法估算的相同。

(我们已经知道 cache-line 拆分存储将需要重放,这会影响吞吐量,因此我们知道即使没有缓存未命中,这也不会完全管理每个结构 1.5 个周期。但这可能是仍然比 Harold 的每 2 个结构 4 个周期 = 每个结构 2 个周期的瓶颈要好。这个速度实际上应该是可以实现的,因为它在不需要在 cache-line 拆分上重播的洗牌上出现瓶颈。)

我预计 Ryzen 的吞吐量会类似,但瓶颈在于存储吞吐量。我们主要使用 128 位向量,Ryzen 的洗牌吞吐量比英特尔更好。在 SnB-family 上,循环中有 4 个 shuffle 微指令。

如果我能以不同的方式洗牌,这样我就可以得到两个连续的结构作为这对向量的高半部分,这将打开将 2 个标量赋值组合成的可能性一个 _mm_storeu_si128 与两个 _mm_storeh_pi (movhps) 64 位存储重叠。 (仍在为其他两个输出结构进行两次混合。)这将使它总共减少到 5 个存储。

但是 shufps 对从何处获取源数据有限制,因此您不能使用它来模拟 unpcklps 或以不同方式交错。

可能值得对 B1 A1 结构使用 palignr,花费额外的 shuffle uop 来保存商店。

我没有对此进行基准测试,也没有计算未对齐的存储跨越高速缓存行边界的频率(因此会影响吞吐量)。


AVX512

如果我们有 AVX512,我们会有 2 输入 lane-crossing 洗牌,可以让我们更有效地构建 float+int 数据的向量,用更少的每个结构随机播放和存储指令。 (我们可以将 vpermt2ps 与 merge-masking 一起使用到 set1(integer) 中,以便在正确的位置交错 2 个转换结果向量和整数。)

Intel's 4x3 transposition example 启发并基于@PeterCordes 解决方案,这是一个 AVX1 解决方案,它应该在 8 个周期内获得 8 个结构的吞吐量(瓶颈仍然是 p5):

#include <immintrin.h>
#include <stddef.h>

struct f2u { 
  float O1, O2;
  unsigned int Offset;
};
static const unsigned uiDefaultOffset = 123;

void cvt_interleave_avx(f2u *__restrict dst, double *__restrict pA, double *__restrict pB, ptrdiff_t len)
{
    __m256 voffset = _mm256_castsi256_ps(_mm256_set1_epi32(uiDefaultOffset));

    // 8 structs per iteration
    ptrdiff_t i=0;
    for(; i<len-7; i+=8)
    {
        // destination address for next 8 structs as float*:
        float* dst_f = reinterpret_cast<float*>(dst + i);

        // 4*vcvtpd2ps    --->  4*(p1,p5,p23)
        __m128 inA3210 = _mm256_cvtpd_ps(_mm256_loadu_pd(&pA[i]));
        __m128 inB3210 = _mm256_cvtpd_ps(_mm256_loadu_pd(&pB[i]));
        __m128 inA7654 = _mm256_cvtpd_ps(_mm256_loadu_pd(&pA[i+4]));
        __m128 inB7654 = _mm256_cvtpd_ps(_mm256_loadu_pd(&pB[i+4]));

        // 2*vinsertf128  --->  2*p5
        __m256 A76543210 = _mm256_set_m128(inA7654,inA3210);
        __m256 B76543210 = _mm256_set_m128(inB7654,inB3210);

        // 2*vpermilps    --->  2*p5
        __m256 A56741230 = _mm256_shuffle_ps(A76543210,A76543210,_MM_SHUFFLE(1,2,3,0));
        __m256 B67452301 = _mm256_shuffle_ps(B76543210,B76543210,_MM_SHUFFLE(2,3,0,1));

        // 6*vblendps     ---> 6*p015 (does not need to use p5)
        __m256 outA1__B0A0 = _mm256_blend_ps(A56741230,B67452301,2+16*2);
        __m256 outA1ccB0A0 = _mm256_blend_ps(outA1__B0A0,voffset,4+16*4);

        __m256 outB2A2__B1 = _mm256_blend_ps(B67452301,A56741230,4+16*4);
        __m256 outB2A2ccB1 = _mm256_blend_ps(outB2A2__B1,voffset,2+16*2);

        __m256 outccB3__cc = _mm256_blend_ps(voffset,B67452301,4+16*4);
        __m256 outccB3A3cc = _mm256_blend_ps(outccB3__cc,A56741230,2+16*2);

        // 3* vmovups     ---> 3*(p237,p4)
        _mm_storeu_ps(dst_f+ 0,_mm256_castps256_ps128(outA1ccB0A0));
        _mm_storeu_ps(dst_f+ 4,_mm256_castps256_ps128(outB2A2ccB1));
        _mm_storeu_ps(dst_f+ 8,_mm256_castps256_ps128(outccB3A3cc));
        // 3*vextractf128 ---> 3*(p23,p4)
        _mm_storeu_ps(dst_f+12,_mm256_extractf128_ps(outA1ccB0A0,1));
        _mm_storeu_ps(dst_f+16,_mm256_extractf128_ps(outB2A2ccB1,1));
        _mm_storeu_ps(dst_f+20,_mm256_extractf128_ps(outccB3A3cc,1));
    }

    // scalar cleanup for  if _iNum is not even
    for (; i < len; i++)
    {
        dst[i].O1 = static_cast<float>(pA[i]);
        dst[i].O2 = static_cast<float>(pB[i]);
        dst[i].Offset = uiDefaultOffset;
    }
}

Godbolt link,最后有最小 test-code:https://godbolt.org/z/0kTO2b

出于某种原因,gcc 不喜欢生成直接从内存转换为寄存器的 vcvtpd2ps。这 可能 works better with aligned loads(无论如何,输入和输出对齐可能是有益的)。 clang 显然想在最后用 vextractf128 指令之一打败我。