使用 simd 在双打数组中找到 nan

Question

这个问题非常类似于：

用于浮点相等比较的 SIMD 指令（使用 NaN == NaN）

尽管该问题侧重于 128 位向量，并且对识别 +0 和 -0 有要求。

我有一种感觉，我自己可能会得到这个，但英特尔内在指南页面似乎已关闭：/

我的目标是获取一个双精度数组并返回数组中是否存在 NaN。我预计大多数情况下不会有一条路线，并希望该路线具有最佳性能。

最初，我打算对自己进行 4 个双打的比较，反映用于 NaN 检测的非 SIMD 方法（即 NaN 只有值为a != atrue 的值）。就像是：

data *double = ...
__m256d a, b;
int temp = 0;

//This bit would be in a loop over the array
//I'd probably put a sentinel in and loop over while !temp
a = _mm256_loadu_pd(data);
b = _mm256_cmp_pd(a, a, _CMP_NEQ_UQ);
temp = temp | _mm256_movemask_pd(b);

但是，在一些比较示例中，除了比较本身之外，似乎还进行了某种 NaN 检测。我简单地想，如果像这样的东西 _CMP_EQ_UQ可以检测到 NaN，我可以使用它，然后我可以比较 4 个双打和 4 个双打，并神奇地同时查看 8 个双打。

__m256d a, b, c;
a = _mm256_loadu_pd(data);
b = _mm256_loadu_pd(data+4);
c = _mm256_cmp_pd(a, b, _CMP_EQ_UQ);

在这一点上，我意识到我并不是很直接，因为我可能碰巧将一个不是 NaN 的数字与其自身进行比较（即 3 == 3）并以这种方式得到命中。

所以我的问题是，将 4 个双精度数与自己进行比较（如上所述）是我能做的最好的事情，还是有其他更好的方法来确定我的数组是否有 NaN？

Answer 1

您可能可以通过检查 fenv 状态来完全避免这种情况，或者如果没有，则缓存阻止它和/或将其折叠到另一个传递相同数据的通道中，因为它的计算强度非常低（加载/存储的每个字节的工作量），所以它容易成为内存带宽的瓶颈。见下文。

您要查找的比较谓词是_CMP_UNORD_Qor_CMP_ORD_Q告诉您比较是无序或有序的，即至少一个操作数是 NaN，或者两个操作数分别是非 NaN。 有序/无序比较是什么意思？

用于cmppd列出谓词的 asm 文档与内在指南具有相同或更好的细节。

所以是的，如果您希望 NaN 很少见并希望快速扫描大量非 NaN 值，您可以使用vcmppd两个不同的向量相互对立。如果您关心 NaN 的位置，那么一旦您知道两个输入向量中的任何一个中至少有一个，您就可以做额外的工作来解决这个问题。（喜欢_mm256_cmp_pd(a,a, _CMP_UNORD_Q)为最低设置位提供 movemask + bitscan。）

OR 或 AND 多次比较每个 movemask

与其他 SSE/AVX 搜索循环一样，您也可以movemask通过将一些比较结果与_mm256_or_pd（查找任何无序）或_mm256_and_pd（检查所有有序）组合来分摊成本。例如，检查每个移动掩码/测试/分支的几个缓存行（4x_mm256d和 2x _mm256_cmp_pd）。 （glibc 的 asmmemchr并strlen使用这个技巧。）同样，这优化了您不希望提前退出并且必须扫描整个阵列的常见情况。

还要记住，检查相同的元素两次完全没问题，所以你的清理可以很简单：一个加载到数组末尾的向量，可能与你已经检查过的元素重叠。

// checks 4 vectors = 16 doubles
// non-zero means there was a NaN somewhere in p[0..15]
static inline
int any_nan_block(double *p) {
    __m256d a = _mm256_loadu_pd(p+0);
    __m256d abnan = _mm256_cmp_pd(a, _mm256_loadu_pd(p+ 4), _CMP_UNORD_Q);
    __m256d c = _mm256_loadu_pd(p+8);
    __m256d cdnan = _mm256_cmp_pd(c, _mm256_loadu_pd(p+12), _CMP_UNORD_Q);
    __m256d abcdnan = _mm256_or_pd(abnan, cdnan);
    return _mm256_movemask_pd(abcdnan);
}
// more aggressive ORing is possible but probably not needed
// especially if you expect any memory bottlenecks.

我把 C 写成汇编，每个源代码行一条指令。（加载/内存源 cmppd）。如果在 Intel 上使用非索引寻址模式，那么这 6 条指令在现代 CPU 的融合域中都是单指令的。 test/jnz作为break条件将其提高到 7 uops。

在循环中，add reg, 16*8指针增量是另一个 1 uop，并且cmp / jne作为循环条件又增加一个，使其达到 9 uop。所以不幸的是，在 Skylake 上，前端的瓶颈为 4 uop/时钟，至少需要 9/4 个周期来发布 1 次迭代，并没有完全饱和加载端口。Zen 2 或 Ice Lake 可以在每个时钟承受 2 个负载，而无需更多展开或其他级别的vorpd组合。

另一个可能的技巧是在两个向量上使用vptest或vtestpd来检查它们是否都非零。但我不确定是否可以正确检查两个向量的每个元素都非零。 是否可以使用 PTEST 来测试两个寄存器是否都为零或其他条件？表明另一种方式（_CMP_UNORD_Q输入都为零）是不可能的。

但这并没有真正的帮助：vtestpd/jcc总共是 3 个 uops，而vorpd/ vmovmskpd/test+jcc在具有 AVX 的现有 Intel/AMD CPU 上也是 3 个融合域 uops，因此当您对结果进行分支时，它甚至不是吞吐量的胜利. 因此，即使有可能，它也可能收支平衡，尽管它可能会节省一些代码大小。如果从全 1 的情况中挑选出全零或 mix_zeros_and_ones 的情况需要多个分支，则不值得考虑。

避免工作：fenv改为检查标志

如果您的数组是此线程中的计算结果，只需检查 FP 异常粘滞标志（在 MXCSR 中手动或通过fenv.h fegetexcept）以查看自您上次清除 FP 异常以来是否发生了 FP“无效”异常。如果没有，我认为这意味着 FPU 没有产生任何 NaN 输出，因此此线程从那时起写入的数组中没有任何输出。

如果已设置，则必须检查；对于未传播到此数组的临时结果，可能已引发无效异常。

缓存阻塞：

如果/当 fenv 标志不能让您完全避免工作，或者不是您程序的好策略，请尝试将此检查折叠到生成数组的任何内容中，或折叠到读取它的下一次传递中。因此，当数据已经加载到向量寄存器中时，您正在重用数据，从而增加计算强度。（每个加载/存储的 ALU 工作。）

即使数据在 L1d 中已经很热，它仍然会在加载端口带宽上出现cmppd瓶颈：在 2/clock 加载端口带宽上，在具有 2/clock 的 CPU vcmppd ymm（Skylake 但不是 Haswell）上，每个仍然存在2 个加载瓶颈。

同样值得对齐您的指针以确保您从 L1d 缓存中获得完整的负载吞吐量，尤其是当 L1d 中的数据有时已经很热时。

或者至少 缓存块它，以便您在缓存中很热时在同一块上运行另一个循环之前检查一个 128kiB 块。这是 256k L2 大小的一半，因此您的数据在上一次传递中应该仍然很热，和/或在下一次传递中仍然很热。

绝对避免在整个多兆字节阵列上运行它并支付将其从 DRAM 或 L3 缓存进入 CPU 内核的成本，然后在另一个循环读取它之前再次驱逐。这是最坏情况下的计算强度，需要付出多次将其放入 CPU 内核的私有缓存的成本。