c++ - OpenMP 的 SIMD 指令可以矢量化索引操作吗？

Question

假设我有一个 MxN 矩阵 (SIG) 和一个 Nx1 分数索引列表 (idxt)。 idxt 中的每个小数索引唯一对应于 SIG 中的相同位置列。我想使用 idxt 中存储的索引来索引 SIG 中的适当值，获取该值并将其保存在另一个 Nx1 vector 中。由于 idxt 中的索引是小数，因此我需要在 SIG 中进行插值。这是使用线性插值的实现:

void calcPoint(const Eigen::Ref<const Eigen::VectorXd>& idxt,
               const Eigen::Ref<const Eigen::Matrix<short int, -1, -1>>& SIG,
               double& bfVal) {
    Eigen::VectorXd bfPTVec(idxt.size());
    #pragma omp simd
    for (int j = 0; j < idxt.size(); j++) {
        int vIDX = static_cast<int>(idxt(j));
        double interp1 = vIDX + 1 - idxt(j);
        double interp2 = idxt(j) - vIDX;
        bfPTVec(j) = (SIG(vIDX,j)*interp1 + SIG(vIDX+1,j)*interp2);
    }
    bfVal = ((idxt.array() > 0.0).select(bfPTVec,0.0)).sum();
}

我怀疑这里有更好的方法来实现循环体，这将有助于编译器更好地利用 SIMD 操作。例如，据我了解，强制编译器在类型之间进行转换(无论是像第一行那样显式转换还是像某些数学运算那样隐式转换)都不是可矢量化操作。

此外，通过使对 SIG 的访问依赖于运行时计算的 idxt 中的值，我不清楚我在这里执行的内存读写类型是否可向量化，或者如何向量化。看看我的问题的大图描述，其中每个 idxt 对应于与 SIG 相同的“位置”列，我感觉到它应该是一个可矢量化的操作，但我不确定如何将其转换为好的代码。

澄清

感谢这些评论，我意识到我没有指定当 idxt 在此方法之外初始化时，我不希望对 idxt 中的最终求和做出贡献的某些值设置为零。因此，上面给出的示例中的最后一行。

Answer 1

理论上，假设处理器支持此操作，它应该是可能的。但实际上，由于多种原因，情况并非如此。

首先，支持指令集 AVX-2(或 AVX-512)的主流 x86-64 处理器确实有相关指令:gather SIMD instructions 。不幸的是，指令集非常有限:您只能从基于 32 位/64 位索引的内存中获取 32 位/64 位值。此外，该指令在主流处理器上还没有非常有效地实现。事实上，它单独获取每个项目，这并不比标量代码快，但如果代码的其余部分被矢量化，这仍然很有用，因为手动读取许多标量值来填充 SIMD 寄存器往往效率较低(尽管由于收集指令的早期实现效率相当低，它在旧处理器上的速度出奇地快)。请注意，SIG 矩阵很大，那么缓存未命中将显着减慢代码速度。

此外，主流处理器上默认不启用 AVX-2，因为并非所有 x86-64 处理器都支持它。因此，您需要启用 AVX-2(例如使用 -mavx2)，以便编译器可以有效地矢量化循环。不幸的是，这还不够。事实上，大多数编译器目前无法自动检测何时可以/应该使用该指令。即使他们可以，IEEE-754 float 运算不具有关联性并且值可以是无穷大或 NaN，这一事实通常不会帮助他们生成有效的代码(尽管这里应该没问题)。请注意，您可以告诉编译器可以假定操作是关联的，并且您仅使用有限/基本实数(例如使用 -ffast-math，这可能是不安全的)。如果编译器无法完全内联所有函数(ICC 就是这种情况)，同样的情况也适用于 Eigen 类型/运算符。

为了加快代码速度，您可以尝试将 SIG 变量的类型更改为包含 int32_t 项的矩阵引用。另一种可能的优化是将循环拆分为固定大小的小块(例如 32 个项目)，并将循环拆分为许多部分，以便在单独的循环中计算间接寻址，以便编译器可以对至少某些循环进行矢量化。一些编译器(例如 Clang)能够自动为您执行此操作:它们为循环的一部分生成快速 SIMD 实现，并使用标量指令进行间接寻址。如果这还不够(到目前为止似乎是这种情况)，那么您当然需要使用 SIMD 内在函数自行对循环进行矢量化(或者可能使用可以为您完成此操作的 SIMD 库)。