代码之家  ›  专栏  ›  技术社区  ›  user276648 Jeremy Cook

在C++ SIMD中转换签名短到浮点

  •  4
  • user276648 Jeremy Cook  · 技术社区  · 6 年前

    我有一个有符号短数组,我想除以2048,得到一个浮点数组。

    我发现 SSE: convert short integer to float 允许转换 未签名 短裤要浮起来,但我也要处理签名短裤。

    下面的代码可以工作,但只适用于正短路。

    // We want to divide some signed short by 2048 and get a float.
    const auto floatScale = _mm256_set1_ps(2048);
    
    short* shortsInput = /* values from somewhere */;
    float* floatsOutput = /* initialized */;
    
    __m128i* m128iInput = (__m128i*)&shortsInput[0];
    
    // Converts the short vectors to 2 float vectors. This works, but only for positive shorts.
    __m128i m128iLow = _mm_unpacklo_epi16(m128iInput[0], _mm_setzero_si128());
    __m128i m128iHigh = _mm_unpackhi_epi16(m128iInput[0], _mm_setzero_si128());
    __m128 m128Low = _mm_cvtepi32_ps(m128iLow);
    __m128 m128High = _mm_cvtepi32_ps(m128iHigh);
    
    // Puts the 2 __m128 vectors into 1 __m256.
    __m256 singleComplete = _mm256_castps128_ps256(m128Low);
    singleComplete = _mm256_insertf128_ps(singleComplete, m128High, 1);
    
    // Finally do the math
    __m256 scaledVect = _mm256_div_ps(singleComplete, floatScale);
    
    // and puts the result where needed.
    _mm256_storeu_ps(floatsOutput[0], scaledVect);
    

    我怎样才能把我签名的短裤换成花边呢?或者也许有更好的方法来解决这个问题?


    编辑: 与非SIMD算法相比,我尝试了不同的答案,在2048阵列上,在AMD Ryzen7 2700上,在3.2GHz的频率下,做了10倍的运算。我使用的是Visual15.7.3,主要是默认配置:

    /permissive- /Yu"stdafx.h" /GS /GL /W3 /Gy /Zc:wchar_t /Zi /Gm- /O2 /sdl 
    /Fd"x64\Release\vc141.pdb" /Zc:inline /fp:precise /D "NDEBUG" /D "_CONSOLE"
    /D "_UNICODE" /D "UNICODE" /errorReport:prompt /WX- /Zc:forScope
    /arch:AVX2 /Gd /Oi /MD /openmp /FC /Fa"x64\Release\" /EHsc /nologo
    /Fo"x64\Release\" /Fp"x64\Release\test.pch" /diagnostics:classic 
    

    请注意,我对SIMD非常陌生,多年来一直没有使用C++。以下是我得到的结果(我分别地重新运行每个测试,而不是一个接一个地重新运行,得到了这样更好的结果):

    • 无simd:7300ms
    • WIM的答案:2300ms
    • CHTZ的SSE2回答:1650毫秒
    • CHTZ的AVX2答案:2100ms

    因此,我通过使用simd得到了一个很好的加速,而chtz的sse2答案虽然更加冗长和复杂,但速度更快。(至少在启用了AVX的情况下编译时,它避免了使用3操作数的VEX编码指令复制寄存器的额外指令。在英特尔CPU上,AVX2版本应该比128位版本快得多。)

    这是我的测试代码:

    const int size = 2048;
    const int loopSize = (int)1e7;
    
    float* noSimd(short* shortsInput) {
        float* floatsOutput = new float[size];
    
        auto startTime = std::chrono::high_resolution_clock::now();
    
        for (int i = 0; i < loopSize; i++) {
            for (int j = 0; j < size; j++) {
                floatsOutput[j] = shortsInput[j] / 2048.0f;
            }
        }
    
        auto stopTime = std::chrono::high_resolution_clock::now();
        long long totalTime = (stopTime - startTime).count();
    
        printf("%lld noSimd\n", totalTime);
    
        return floatsOutput;
    }
    
    float* wimMethod(short* shortsInput) {
        const auto floatScale = _mm256_set1_ps(1.0f / 2048.0f);
        float* floatsOutput = new float[size];
    
        auto startTime = std::chrono::high_resolution_clock::now();
    
        for (int i = 0; i < loopSize; i++) {
            for (int j = 0; j < size; j += 8) {
                __m128i short_vec = _mm_loadu_si128((__m128i*)&shortsInput[j]);
                __m256i int_vec = _mm256_cvtepi16_epi32(short_vec);
                __m256  singleComplete = _mm256_cvtepi32_ps(int_vec);
    
                // Finally do the math
                __m256 scaledVect = _mm256_mul_ps(singleComplete, floatScale);
    
                // and puts the result where needed.
                _mm256_storeu_ps(&floatsOutput[j], scaledVect);
            }
        }
    
        auto stopTime = std::chrono::high_resolution_clock::now();
        long long totalTime = (stopTime - startTime).count();
    
        printf("%lld wimMethod\n", totalTime);
    
        return floatsOutput;
    }
    
    float* chtzMethodSSE2(short* shortsInput) {
        float* floatsOutput = new float[size];
    
        auto startTime = std::chrono::high_resolution_clock::now();
    
        for (int i = 0; i < loopSize; i++) {
            for (int j = 0; j < size; j += 8) {
                // get input:
                __m128i val = _mm_loadu_si128((__m128i*)&shortsInput[j]);
                // add 0x8000 to wrap to unsigned short domain:
                val = _mm_add_epi16(val, const0x8000);
                // interleave with upper part of float(1<<23)/2048.f:
                __m128i lo = _mm_unpacklo_epi16(val, const0x4580);
                __m128i hi = _mm_unpackhi_epi16(val, const0x4580);
                // interpret as float and subtract float((1<<23) + (0x8000))/2048.f
                __m128 lo_f = _mm_sub_ps(_mm_castsi128_ps(lo), constFloat);
                __m128 hi_f = _mm_sub_ps(_mm_castsi128_ps(hi), constFloat);
                // store:
                _mm_storeu_ps(&floatsOutput[j], lo_f);
                _mm_storeu_ps(&floatsOutput[j] + 4, hi_f);
            }
        }
    
        auto stopTime = std::chrono::high_resolution_clock::now();
        long long totalTime = (stopTime - startTime).count();
    
        printf("%lld chtzMethod\n", totalTime);
    
        return floatsOutput;
    }
    
    float* chtzMethodAVX2(short* shortsInput) {
        const auto floatScale = _mm256_set1_ps(1.0f / 2048.0f);
        float* floatsOutput = new float[size];
    
        auto startTime = std::chrono::high_resolution_clock::now();
    
        for (int i = 0; i < loopSize; i++) {
            for (int j = 0; j < size; j += 8) {
    
                // get input:
                __m128i val = _mm_loadu_si128((__m128i*)&shortsInput[j]);
                // interleave with 0x0000
                __m256i val_unpacked = _mm256_cvtepu16_epi32(val);
    
                // 0x4580'8000
                const __m256 magic = _mm256_set1_ps(float((1 << 23) + (1 << 15)) / 2048.f);
                const __m256i magic_i = _mm256_castps_si256(magic);
    
                /// convert by xor-ing and subtracting magic value:
                // VPXOR avoids port5 bottlenecks on Intel CPUs before SKL
                __m256 val_f = _mm256_castsi256_ps(_mm256_xor_si256(val_unpacked, magic_i));
                __m256 converted = _mm256_sub_ps(val_f, magic);
                // store:
                _mm256_storeu_ps(&floatsOutput[j], converted);
            }
        }
    
        auto stopTime = std::chrono::high_resolution_clock::now();
        long long totalTime = (stopTime - startTime).count();
    
        printf("%lld chtzMethod2\n", totalTime);
    
        return floatsOutput;
    }
    
    2 回复  |  直到 6 年前
        1
  •  6
  •   chtz    6 年前

    您可以替换转换epi16->epi32->浮点和乘的标准方法 1.f/2048.f ,通过手动组合浮点。

    这是因为除数是2的幂,所以手动组合浮点表示不同的指数。

    多亏了@petercordes,这里有一个优化的avx2版本,使用xor在翻转整数值的符号位的同时设置32位浮点的高位字节。fp sub将尾数的低位转换为正确的fp值:

    // get input:
    __m128i val = _mm_loadu_si128((__m128i*)input);
    // interleave with 0x0000
    __m256i val_unpacked = _mm256_cvtepu16_epi32(val);
    
    // 0x4580'8000
    const __m256 magic = _mm256_set1_ps(float((1<<23) + (1<<15))/2048.f);
    const __m256i magic_i = _mm256_castps_si256(magic);
    
    /// convert by xor-ing and subtracting magic value:
    // VPXOR avoids port5 bottlenecks on Intel CPUs before SKL
    __m256 val_f = _mm256_castsi256_ps(_mm256_xor_si256(val_unpacked, magic_i));
    __m256 converted = _mm256_sub_ps(val_f, magic);
    // store:
    _mm256_storeu_ps(output, converted);
    

    看到它 on the Godbolt compiler explorer with gcc and clang ;在Skylake i7-6700K上,一个2048元素的高速缓存循环需要大约360个时钟周期,与执行标准符号扩展/转换/乘法(循环展开量相似)的@wim版本的速度(到测量误差内)相同。由@petercordes在linux下测试 perf . 但在Ryzen,这可能会更快,因为我们避免 _mm256_cvtepi32_ps (Ryzen每2个时钟有1个吞吐量 vcvtdq2ps ymm : http://agner.org/optimize/ )

    的异或 0x8000 下半部分等于加/减 0x8000 ,因为忽略溢出/进位。巧合的是,这允许使用相同的魔法常数进行异或运算和减法运算。

    奇怪的是,gcc和clang更喜欢用 -magic ,它不会重复使用常量…他们更喜欢使用 add 因为它是可交换的,但在这种情况下没有好处,因为它们没有将它与内存操作数一起使用。


    这里有一个SSE2版本,它分别执行有符号/无符号翻转,而不设置32位FP位模式的高2字节。

    我们在用一个 _mm_add_epi16 _mm_unpackXX_epi16 还有两个 _mm_sub_ps 对于8个值 _mm_castsi128_ps 没有行动,而且 _mm_set 将缓存在寄存器中):

    // get input:
    __m128i val = _mm_loadu_si128((__m128i*)input);
    // add 0x8000 to wrap to unsigned short domain:
    // val = _mm_add_epi16(val, _mm_set1_epi16(0x8000));
    val = _mm_xor_si128(val, _mm_set1_epi16(0x8000));  // PXOR runs on more ports, avoids competing with FP add/sub or unpack on Sandybridge/Haswell.
    
    // interleave with upper part of float(1<<23)/2048.f:
    __m128i lo = _mm_unpacklo_epi16(val, _mm_set1_epi16(0x4580));
    __m128i hi = _mm_unpackhi_epi16(val, _mm_set1_epi16(0x4580));
    // interpret as float and subtract float((1<<23) + (0x8000))/2048.f
    __m128 lo_f = _mm_sub_ps(_mm_castsi128_ps(lo), _mm_set_ps1(float((1<<23) + (1<<15))/2048.f));
    __m128 hi_f = _mm_sub_ps(_mm_castsi128_ps(hi), _mm_set_ps1(float((1<<23) + (1<<15))/2048.f));
    // store:
    _mm_storeu_ps(output, lo_f);
    _mm_storeu_ps(output+4, hi_f);
    

    使用演示: https://ideone.com/b8BfJd

    如果您的输入是无符号短 , the _增加16毫米 不需要(以及 1<<15 _毫米潜艇 当然需要移除)。那你就有马拉的答案了 SSE: convert short integer to float .

    这个可以 容易地 移植到avx2,每次迭代的转换次数是avx2的两倍,但必须注意输出元素的顺序(感谢@wim指出这一点)。


    此外,对于纯SSE解决方案,可以简单地使用 _mm_cvtpi16_ps ,但这是一个Intel库函数。没有一条指令能做到这一点。

    // cast input pointer:
    __m64* input64 = (__m64*)input;
    // convert and scale:
    __m128 lo_f = _mm_mul_ps(_mm_cvtpi16_ps(input64[0]), _mm_set_ps1(1.f/2048.f));
    __m128 hi_f = _mm_mul_ps(_mm_cvtpi16_ps(input64[1]), _mm_set_ps1(1.f/2048.f));
    

    我没有对任何解决方案进行基准测试(也没有检查理论吞吐量或延迟)

        2
  •  5
  •   Peter Cordes    6 年前

    对于avx2,不需要分别转换高、低部件:

    const auto floatScale = _mm256_set1_ps(1.0f/2048.0f);
    
    short* shortsInput = /* values from somewhere */;
    float* floatsOutput = /* initialized */;
    
    __m128i short_vec = _mm_loadu_si128((__m128i*)shortsInput);
    __m256i int_vec =  _mm256_cvtepi16_epi32 (short_vec);
    __m256  singleComplete = _mm256_cvtepi32_ps (int_vec);
    
    // Finally do the math
    __m256 scaledVect = _mm256_mul_ps(singleComplete, floatScale);
    
    // and puts the result where needed.
    _mm256_storeu_ps(floatsOutput, scaledVect);
    

    这个编译得很好 on the Godbolt compiler explorer ,使用L1D缓存中的热输入/输出和对齐的输入/输出数组,在Skylake i7-6700k上以约360个时钟周期转换2048个元素的数组(在重复循环中测试)。这是每个元素大约0.18个周期,或者每个时钟周期大约5.7个转换。或每个向量约1.4个周期,包括存储。它在前端吞吐量(每时钟3.75个融合域UOP)上大多是瓶颈,甚至在Clang的循环展开时也是如此,因为转换是5 UOP。

    注意 vpmovsxwd ymm, [mem] 即使在haswell/skylake上使用简单的寻址模式,也不能将微融合为单个UOP,因此在这种情况下,最近的gcc/clang转换指针会增加到使用单个循环计数器的索引寻址,这是好事。使用大多数内存源矢量指令(如 vpmovsxwd xmm, [mem] ,这将需要额外的UOP: Micro fusion and addressing modes

    有了一个负载和一个存储,这些存储就不能在haswell/skylake的port7存储agu上运行了,后者只处理非索引寻址模式。

    对于英特尔CPU上的最大吞吐量(如果没有内存瓶颈),需要展开循环,因为加载+转换+存储已经是4个Uops。与@chtz的答案相同。

    理想情况下,如果只需要读取浮点值几次,就可以立即使用向量结果进行进一步的计算。它只有3条指令(但对于无序的执行隐藏有一定的延迟)。在需要时重做转换可能比使用更大的缓存占用空间来存储两倍大的 float[] 产生内存;这取决于您的用例和硬件。