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

Dr.*_*ABT 7 c++ x86 simd intrinsics avx

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

是否有使用SIMD Intrinsics进行复制和转换操作的更快方法?我已经在更快的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;
        }
Run Code Online (Sandbox Code Playgroud)

注意:结构格式为[Float,Float,Int](24字节),但是我们可以(如果有帮助的话)添加额外的4字节填充,使其成为32字节。

har*_*old 7

这是使用SSE4.1的尝试,没有AVX(这样做比较棘手,到目前为止,我提出了更多的改组方法),并使用12byte / point格式:(未测试)

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;
    }
}
Run Code Online (Sandbox Code Playgroud)

这利用了shufps从两个不同的源中进行选择以进行动态数据和恒定偏移量合并的能力,相同的混洗也会移动需要移动的每个组中的浮点。然后使用混合将单个浮动替换为已经在正确位置的另一个浮动。这需要2次混洗和2次混和,还有3次混洗和0次混和的方法,但是混洗都将在当前的Intel处理器上分配到p5,而混和可以转到其他端口。转换也已经使用了p5,因此会陷入困境,使用混合应该会更好。每次迭代仍为4 p5 µops,因此每个处理的项目至少需要2个周期,这并不好。

主循环会跳过最后一项,以便不会超出范围,它会略微重叠16个字节的存储,并在结构末尾写入4个字节。该部分将被下一个存储区的实际结果覆盖,但是在数组末尾执行该操作可能很危险。

  • 256位vcvtpd2ps xmm0和ymmword [mem]的AVX可能在这里很有用,可以将[A3 A2 A1 A0]和B转换为__m128向量,您可以用2种不同的方式一起洗牌(unpcklps和`unpckhps`)。这将为FP零件的`vmovlps`和`vmovhps`存储设置。`vmovhps [mem],xmm`不需要花费Intel SnB系列上的洗牌,只是一家纯粹的商店。但是您确实需要分别复制整数元素,可能使用标量加载/存储来避免随机端口瓶颈。瓶颈=每个输出结构2次存储(直到Ice Lake) (2认同)

Pet*_*des 6

此问题与memcpy不太相似。这一切都是为了优化循环不变整数成员的混洗和/或标量存储的交织。这使SIMD变得很难。

您是否需要这种intfloat成员交错的存储格式?交叉浮动是很糟糕的。我假定以后的代码将以int不同的结构修改s,否则对于每个元素都将其复制是没有意义的。

您是否可以按4个元素的组进行工作,例如 struct { float a[4], b[4]; int i[4]; };可以将4x连续的加载和转换double为4x float并进行128位SIMD存储?访问单个输出数组“结构”的所有3个成员时,您仍然会有一些空间上的局限性。


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

否则,另一种策略是对浮点数与分别使用单独的存储int,因此您不需要重叠。我们可能可以对其进行优化,以至于每个时钟周期1个结构每2个周期的瓶颈将成为瓶颈。(或稍低一些,因为IIRC高速缓存拆分存储至少在Intel CPU上需要重播存储uop。)我们还可以按4*12 = 3*16字节展开并使用SIMD存储来保存2个整数存储,该存储与浮点数据重叠。48个字节= 作为四个结构的一部分xyIx|yIxy|IxyI包含四个I元素,但是它们足够接近,因此我们可以使用两个_mm_storeu_si128( set1(offset) )内在函数存储所有四个元素。然后存储xy与之重叠的对。16字节边界标有|。如果缓存线拆分是一个问题,我们可以做2个标量和其中最后一个向量SIMD 对齐(如果输出阵列排列的16个字节)。或者在Intel Haswell和更高版本的CPU上,一个32字节对齐的存储可能会很好。


如果我们不小心的话,我们很容易在Intel CPU上产生瓶颈,特别是在FP哈希只能在端口5上运行的Sandybridge家族(通过Skylake / Coffee Lake的SnB)上,这是瓶颈。这就是为什么我考虑不进行重组每一个结构总共存储1个存储区。

SIMD double-> float转换的成本为2 uop:shuffle + FP-math,因为float是宽度的一半,并且指令将float打包到向量寄存器的底部。

在此,AVX可用于将4 doubles 转换为4 s的SIMD向量float

除此之外,我同意@harold的观点,认为128位向量可能是一个不错的选择。甚至AVX2也没有很好的2输入交叉穿越功能,并且AVX1的功能非常有限。因此,我们可以使用256位-> 128位double-> float转换来提供基于的交错策略__m128

vmovhps [mem], xmm无需在Intel CPU上进行洗牌,只是一个纯粹的存储,因此将2个向量洗牌在一起,并[ B1 A1 B0 A0 ]进入一个向量,就为我们提供了两个低位和高位的64位存储,而没有任何额外的位洗。

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


或者这是另一个使用上述方法但又进行一些混合以保存存储的想法

在编辑@harold的代码以实现我在上面的文本中所写的想法时,我基本上想到了这一点。在此处使用混合是减少存储和随机端口压力的好方法。

上面的一些想法仍然值得探索,尤其是进行大量存储,set1(offset)然后与64位vmovlps存储重叠。(将3x2 = 6或3x4 = 12的输出结构展开后,使其一次转换为4个double的倍数。) 12 * 12 = 144字节,它是16的倍数,而不是32或64,所以我们至少可以知道在哪里我们一直都相对于16字节的边界,但是除非我们展开更多,否则不缓存行。(可能会留下更多需要清理的工作,并使代码大小膨胀。)

#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;
    }
}
Run Code Online (Sandbox Code Playgroud)

Godbolt -O3 -march=skylake上的gcc9.1将主循环编译为前端的19个融合域uops。(这两种vcvtpd2ps指令都无法进行微融合,因为GCC并没有做任何聪明的事情,例如pB相对于寻址,pA以避免对其中一个进行索引寻址。因此,它们各有3个微字:load + convert + shuffle)

但是无论如何,它都会在后端存储上造成瓶颈,即使从4个宽的前端发出每个迭代要花费整个5个周期。

每次迭代有6个存储区(每4个结构),这将使其瓶颈最多达到每6个周期1次迭代,这是存储数据端口/执行单元上的瓶颈。(直到Ice Lake每个时钟可以存储2个存储。)因此,在理论上最好的情况下,每1.5个周期可以实现1个结构,这与我之前对重叠存储的想法进行的估算相同。

(我们已经知道,缓存行拆分存储将需要重播,这会消耗吞吐量,因此,即使没有缓存未命中,我们也知道这并不能很好地管理每个结构1.5个周期。但是,它可能仍然比Harold的瓶颈4更好。每2个结构的周期=每个结构2个周期。尽管如此,该速度实际上应该是可以实现的,因为它在洗牌时会产生瓶颈,不需要在高速缓存行拆分中重播。

我预计Ryzen的吞吐量将类似,成为商店吞吐量的瓶颈。我们主要使用128位向量,并且Ryzen的洗牌吞吐量比Intel高。在SnB系列上,循环中有4个随机码。

如果我可以进行不同的改组,以便可以得到两个连续的结构作为向量对的上半部分,那将有可能将2个标量分配组合为一个_mm_storeu_si128与两个_mm_storeh_pimovhps)64位存储区重叠的结构。(仍然对其他两个输出结构进行两次混合。)这将使总数减少到5。

但是shufps对从何处获取源数据有限制,因此您不能使用它来unpcklps进行不同的仿真或交错。

可能值得palignr在B1 A1结构中使用它,多花一些乱码来节省商店。

我还没有对此进行基准测试,也没有计算未对齐的存储多长时间一次越过缓存行边界(因此会提高吞吐量)。


AVX512

如果我们有AVX512,我们将有2输入道交叉改组,这可以让我们更有效地构建float + int数据向量,并且每个结构的改组和存储指令更少。(我们可以将vpermt2ps合并屏蔽set1(integer)与一起使用,以在正确的位置插入2个转换结果向量和整数。)


cht*_*htz 5

Loosely inspired by Intel's 4x3 transposition example and based on @PeterCordes solution, here is an AVX1 solution, which should get a throughput of 8 structs within 8 cycles (bottleneck is still 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;
    }
}
Run Code Online (Sandbox Code Playgroud)

Godbolt link, with minimal test-code at the end: https://godbolt.org/z/0kTO2b

For some reason, gcc does not like to generate vcvtpd2ps which directly convert from memory to a register. This might works better with aligned loads (having input and output aligned is likely beneficial anyway). And clang apparently wants to outsmart me with one of the vextractf128 instructions at the end.

  • ^^ 这是在我的 PC 上进行测试的,在处理 10 亿个元素时,相同循环的时钟速度为 550 毫秒,而 1300 毫秒。 (2认同)