GCC 不断抱怨 AVX512 函数 _mm512_cvt_roundpd_epi64 的“错误：不正确的舍入操作数”

Posted 2023-02-16

【中文标题】GCC 不断抱怨 AVX512 函数 _mm512_cvt_roundpd_epi64 的“错误：不正确的舍入操作数”

【英文标题】：GCC keeps complaining "error: incorrect rounding operand" for a AVX512 functions _mm512_cvt_roundpd_epi64

【发布时间】：2018-08-24 04:23:52

【问题描述】：

我正在使用 _mm512_cvt_roundpd_epi64 并不断收到编译器错误：

/dump/1/alicpp2/built/gcc-7.3.0-7u2/gcc-7.3.0/lib/gcc/x86_64-pc-linux-gnu/7.3.0/include/avx512dqintrin.h:1574： 14: 错误：不正确的舍入操作数 __R);

这是我的代码：

    #include <iostream>
    #include <immintrin.h>

    void Date64Align(int64_t* dst, int64_t* src, size_t length) 
      constexpr int dop = 512 / 64;
      int64_t starting_epoch_milliseconds_ = 1513728000;
      int32_t granularity_milliseconds_ = 3600;

      __m512i start = _mm512_set1_epi64(starting_epoch_milliseconds_);
      __m512i granularity = _mm512_set1_epi64(granularity_milliseconds_);

      double temp = (double)granularity_milliseconds_;
      __m512d granularity_double = _mm512_set1_pd(temp);

      for (int i = 0; i < length / dop; ++i) 
        // load the src (load X into SIMD register
        __m512i data = _mm512_load_epi64(src);
        // X - starting_epoch_milliseconds_
        data = _mm512_sub_epi64(data, start);
        // convert X to double
        __m512d double_data;
        double_data = _mm512_cvt_roundepi64_pd(data, _MM_FROUND_TO_NEAREST_INT);

        // X = X / Y in double
        double_data = _mm512_div_pd(double_data, granularity_double);

        // Convert X to int64
        data = _mm512_cvt_roundpd_epi64(double_data, _MM_FROUND_NO_EXC);

        data = _mm512_mullo_epi64(data, granularity);

        // store X
        _mm512_store_epi64(dst, data);

        src += dop;
        dst += dop;
      
    

    int main() 
      return 0;
    

还有我的 CMakeFileLists.txt：

    cmake_minimum_required(VERSION 3.11)
    project(untitled3)

    set(CMAKE_CXX_STANDARD 17)
    set(CMAKE_CXX_FLAGS "$CMAKE_CXX_FLAGS -ggdb -msse4.2 -mavx512f - 
    mavx512dq")


    add_executable(untitled3 main.cpp)

有人熟悉 AVX512 库并帮助回答我的问题吗？

【问题讨论】：

你能发一个minimal reproducible example，我实际上可以复制/粘贴到gcc.godbolt.org看看吗？你#include <immintrin.h> 并使用-march=native 或-march=skylake-avx512 编译？

你确定你不能使用 2 的幂粒度，所以你可以移动而不是使用相对较慢的 mullo_epi64 吗？ （并乘以逆而不是除。）不幸的是_mm512_roundscale_pd 不起作用，它不允许负数的小数位。 (instrinsic _mm512_round_ps is missing for AVX512)

感谢您的回复。我用一个最小的例子更新了代码，并附上了一个 CMakeFileLists.txt。请你帮忙看看好吗？

【参考方案1】：

仅供参考，您通常不需要显式四舍五入。默认模式是四舍五入，所有异常都被屏蔽。普通的_mm512_cvtepi64_pd 和_mm512_cvtpd_epi64 的行为与您正在执行的操作相同，除非您使用fenv 或_MM_SET_ROUNDING_MODE 更改了此线程中的默认舍入模式或异常掩码。

抑制异常仅意味着它们没有故障，但如果我正确阅读英特尔的手册，它不会阻止在 MXCSR 中设置相关的粘滞状态位的次正常或溢出。他们说这就像在 MXCSR 中设置了掩码位，并不是说它根本不会阻止在 MXCSR 状态位中记录异常。

_mm512_cvt_roundpd_epi64 更常见的用例是使用 floor 或 ceil 舍入（朝向 -/+Infinity）转换为整数，而不是在转换之前使用单独的舍入步骤 128-位或 256 位向量。

但是，如果您在一些未屏蔽的 FP 异常或可能的非默认舍入模式下运行，那么显式舍入到最近值确实有意义。

---

舍入模式覆盖必须始终包含_MM_FROUND_NO_EXC

如果编译器提供更好的错误消息告诉你这一点，那就太好了。 （TODO：提交关于 gcc 和 clang 的功能请求错误报告）。

（_MM_FROUND_CUR_DIRECTION 不计算在内，它的意思是“不覆盖”，就像您使用普通的非 round 版本的内在函数一样。）

英特尔的内在函数指南指出了这一点（在 the entry for _mm512_cvt_roundepi64_pd specifically 中，但您会在每个采用舍入模式覆盖 arg 的内在函数中找到相同的内容。）

根据rounding参数进行舍入，可以是1 的：

(_MM_FROUND_TO_NEAREST_INT |_MM_FROUND_NO_EXC) // round to nearest, and suppress exceptions
(_MM_FROUND_TO_NEG_INF |_MM_FROUND_NO_EXC)     // round down, and suppress exceptions
(_MM_FROUND_TO_POS_INF |_MM_FROUND_NO_EXC)     // round up, and suppress exceptions
(_MM_FROUND_TO_ZERO |_MM_FROUND_NO_EXC)        // truncate, and suppress exceptions
_MM_FROUND_CUR_DIRECTION // use MXCSR.RC; see _MM_SET_ROUNDING_MODE

请注意，_MM_FROUND_NO_EXC 本身恰好是有效的，因为 _MM_FROUND_TO_NEAREST_INT 恰好是 0，这与设置 EVEX.b 时 2 位舍入模式字段的机器代码编码相同。但是你真的应该在_mm512_cvt_roundpd_epi64 中明确说明这一点。

对于没有舍入控制的指令，例如_mm512_cvtt_roundpd_epi64（注意额外的t 用于截断），只允许_MM_FROUND_NO_EXC（或_MM_FROUND_CUR_DIRECTION），因为行为不受2 位字段的值，无论是否指定了舍入覆盖。

---

在 EVEX 前缀的机器编码中，设置舍入模式覆盖意味着 SAE（抑制所有异常）。如果不抑制异常，就无法对 _MM_FROUND_TO_NEAREST_INT 覆盖进行编码。

From Intel's vol.2 instruction set reference manual:

2.6.8 EVEX 中的静态舍入支持

嵌入在 EVEX 编码系统中的静态舍入控制适用 仅用于寄存器到寄存器的浮点指令风格 在两个不同的向量长度处具有舍入语义：（i）标量， (ii) 512 位。在这两种情况下，字段EVEX.L’L 表示舍入 如果设置了EVEX.b，则模式控制覆盖MXCSR.RC。 当设置EVEX.b 时， 隐含“禁止所有异常”。

请注意，舍入覆盖使编译器无法使用内存源操作数，因为在该上下文中，EVEX.b 位表示广播与非广播。

在你的情况下不是问题；数据来自_mm512_sub_epi64，但值得指出的是，在某些情况下不需要额外的加载指令，对已经默认的舍入模式的覆盖可能会产生轻微的性能损失.不过，静态舍入总是比额外的_mm512_roundscale_pd 好（instrinsic _mm512_round_ps is missing for AVX512）。

---

顺便说一句，这些限制（仅适用于标量或 512 位向量，并且仅适用于非内存指令）是 AVX512 有 vcvttpd2qq 的原因，而不是仅仅使用 _MM_FROUND_TO_ZERO|_MM_FROUND_NO_EXC 来表示 _mm512_cvt_roundpd_epi64 .因为没有_mm256_cvt_roundpd_epi64，如果编译器可以将负载折叠到vcvttpd2qq 的内存操作数中，这偶尔会很好。

还有历史先例：自 SSE1 cvttss2si 和 cvttps2dq 以来，英特尔已经进行了截断转换，这使得在不更改 MXCSR 舍入的情况下实现 C 的 FP->int 转换语义更加更加高效以我们过去使用 x87 的方式（在 SSE3 fisttp 之前）。

在 AVX512 之前，从未支持涉及 64 位整数的打包转换，因此该指令不存在现有的 128 位或 256 位版本。不过，提供一个是一个很好的设计决策。

舍入覆盖是 AVX512 中的新功能。在此之前，SSE4.1 roundps / roundpd 可以使用显式模式打包舍入到整数（输入和输出均为 __m128 或 __m128d）。

---

提高效率的替代实现：

添加而不是子：

__m512i minus_start = _mm512_set1_epi64(-starting_epoch_milliseconds_);

 for() 
    __m512i data = _mm512_add_epi64(data, minus_start);
 

add 是可交换的，因此编译器可以将加载折叠成像vpaddq zmm0, zmm8, [rdi] 这样的加载+添加指令，而不是单独的加载+子。 clang 为你做了这个优化，但是gcc doesn't

---

您似乎想将输入整数四舍五入到最接近的 3600 倍数。

用乘法代替除法

1.0/3600 四舍五入到最接近的double 大约是2.777777777777777775368439616699e-04，在 2^53 (the significand precision of double) 中最多只有 0.5 个部分是错误的。大约是 10^-16。对于小于该值的输入，lrint(x * (1.0/3600)) 在lrint(x / 3600.0) 的 1 以内。对于大多数合理大小的输入，它们完全相等。

在乘法之后，你仍然会得到 3600 的 一个 精确倍数，但是在“除法”中有一个微小的错误，你最终可能会偏离 3600 的一倍。

您可以编写一个测试程序来查找除法与乘法运算得到不同结果的情况。

---

您可以将此作为另一次数据传递的一部分吗？对于所有内存带宽来说，计算量并不大。或者，如果您不能将 div_pd 替换为乘以倒数，那么它完全会成为 FP 除法的瓶颈，而不会让其他执行单元保持忙碌。

这里有三种策略：

纯整数，使用乘法逆进行精确除法。 Why does GCC use multiplication by a strange number in implementing integer division?。 Evan AVX512DQ 没有整数乘法，它可以为您提供 64x64 => 128 的 高 一半，只有 vpmullq 64x64 => 64 位（而且它是多个微指令）。

没有 AVX512IFMA VPMADD52HUQ（52x52=>52 位乘法的高半部分），请参阅 Can I use the AVX FMA units to do bit-exact 52 bit integer multiplications?。 （或者，如果您实际上只关心输入的低 32 位，那么 32x32=>64 位乘法和 64 位移位应该可以使用，使用 _mm512_mul_epu32，单微指令 vpmuludq。）但这也需要额外的工作舍入到最近而不是截断。

您现在在做什么：double 除（或乘以倒数），转换为最接近的int64_t，64 位整数相乘。

如果 > 2^53，输入可能会四舍五入到最接近的 double，但最终结果将始终是 3600 的精确倍数（除非乘法溢出 int64_t）。

double 除（或乘），舍入到最接近的整数（不转换），double 乘，转换为整数。

如果最后一次乘法的结果是above 2^(53+4)，则可能会出现问题。 3600 是 2^4 的倍数，但不是 2^5 的倍数。因此，对于非常大的输入，四舍五入到最接近的可表示 double 可能给出的数字不是 3600 的精确倍数。

如果范围限制不是问题，您甚至可以使用fma(val, 3600, -3600.0*start) 将减法折叠起来。

SIMD FP 乘法的吞吐量明显优于整数乘法，因此总体而言它可能是一个胜利，即使有 FP 舍入到最近指令的额外成本。

您有时可以通过添加然后减去一个大常数来避免显式舍入指令，就像 @Mysticial 在 Can I use the AVX FMA units to do bit-exact 52 bit integer multiplications? 中所做的那样。您使值足够大，使其最接近的可表示doubles 是整数。 （How to efficiently perform double/int64 conversions with SSE/AVX?，对于有限范围的输入，也展示了一些 FP 操作技巧。）

也许我们可以rounded=fma(v, 1.0/3600, round_constant)，然后减去round_constant，得到一个不带_mm512_roundscale_pd的四舍五入到最接近整数的值。我们甚至可以使用fma(rounded, 3600, -3600*round_constant) 将其折叠为向上扩展：2^52 * 3600 = 4503599627370496.0 * 3600 完全可以表示为double。

可能存在双舍入问题：首先从int64_t 转换为最接近的double（如果它太大以至于整数不能精确表示），然后在除法和四舍五入到最接近的整数时再次.

---

成本：我假设您可以将 FP 除法替换为乘以 1.0/3600。

1234563 ) = 2 个 FMA 端口的 6 微指令。 vpsubq zmm 也竞争相同的端口，所以真的是 7。 SKX uop 计数来自Agner Fog's testing。

fp 一切：vcvtqq2pd（FMA 端口为 1 uop）+vmulpd（1 uop）+vrndscalepd（2 uop）+vmulpd（1 uop）+vcvtpd2qq（1 uop ) = 6 微秒，但延迟可能更低。 （vrndscale+vmulpd 是 8+4 延迟，比 vpmullq 15 周期延迟快）。但是，如果在数组上循环独立向量，OoO exec 应该很容易隐藏这种延迟，因此节省延迟并不是什么大问题。

我不确定你可以如何高效地进行“整数”乘法，或者使用 FP bithacks 来避免转换指令。如果这对性能至关重要，那可能值得研究。

【讨论】：

非常感谢您的帮助，彼得！我改变了模式，它的工作原理！

@user10267770：顺便说一句，我更新了一些关于提高效率的建议。到目前为止，最重要的是用 mul 替换 div，如果这足够准确的话。因为每 16c 有一个 vdivpd zmm，所以大约是 5 倍 (16/3)。

感谢您进一步的高级帮助，彼得！我用 mul 替换了 div，但这里没有性能提升。甚至还有一点点损失。顺便说一句，我没有将数字四舍五入到最接近的 3600 倍数。我将数字向下舍入到 3600 的倍数（可以是任何其他数字）。这个数字是以毫秒为单位的纪元时间，所以 2^52 对我的一生来说已经足够了。

@user10267770：如果你想向下取整，为什么要使用_mm512_cvt_roundpd_epi64 和_MM_FROUND_NO_EXC（取最近舍入）？我猜您的意思是使用_mm512_cvtt_roundpd_epi64（截断）或_MM_FROUND_TO_ZERO 舍入覆盖？无论如何，您确定您在编译时启用了优化吗？在循环之外计算 1./3600 然后只在内部相乘比你现在做的慢是不合理的。