OpenMP 的 SIMD 指令可以向量化索引操作吗？

Question

假设我有一个 MxN 矩阵 (SIG) 和一个 Nx1 分数索引列表 (idxt)。 idxt 中的每个分数索引唯一对应于 SIG 中的相同位置列。我想使用存储在 idxt 中的索引索引到 SIG 中的适当值，获取该值并将其保存在另一个 Nx1 向量中。由于 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-bit/64-bit 索引的内存中获取 32-bit/64-bit 个值。此外，这条指令在主流处理器上的实现还不是很有效。事实上，它单独获取每个项目并不比标量代码快，但如果代码的其余部分是矢量化的，这仍然有用，因为手动读取许多标量值以填充 SIMD 寄存器的效率往往较低（尽管由于收集指令的早期实现效率非常低，因此在旧处理器上它的速度出奇地快）。请注意，SIG 矩阵很大，那么 缓存未命中 会显着降低代码速度。

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

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

Answer 2

可能不会，但我希望手动矢量化版本更快。

下面是该内部循环的示例，未经测试。它不使用 AVX only SSE up to 4.1，并且应该与你那里的这些特征矩阵兼容。

pIndex 输入指针应指向 idxt 向量的第 j 个元素，pSignsColumn 应指向 SIG 第 j 列的开始矩阵。它假定您的 SIG 矩阵是主要列。它通常是 Eigen 中的默认内存布局，但可以用模板参数覆盖，也可能用宏覆盖。

inline double computePoint( const double* pIndex, const int16_t* pSignsColumn )
{
    // Load the index value into both lanes of the vector
    __m128d idx = _mm_loaddup_pd( pIndex );

    // Convert into int32 with truncation; this assumes the number there ain't negative.
    const int iFloor = _mm_cvttsd_si32( idx );

    // Compute fractional part
    idx = _mm_sub_pd( idx, _mm_floor_pd( idx ) );
    // Compute interpolation coefficients, they are [ 1.0 - idx, idx ]
    idx = _mm_addsub_pd( _mm_set_sd( 1.0 ), idx );

    // Load two int16_t values from sequential addresses
    const __m128i signsInt = _mm_loadu_si32( pSignsColumn + iFloor );
    // Upcast them to int32, then to fp64
    const __m128d signs = _mm_cvtepi32_pd( _mm_cvtepi16_epi32( signsInt ) );

    // Compute the result
    __m128d res = _mm_mul_pd( idx, signs );
    res = _mm_add_sd( res, _mm_unpackhi_pd( res, res ) );
    // The above 2 lines (3 instructions) can be replaced with the following one:
    // const __m128d res = _mm_dp_pd( idx, signs, 0b110001 );
    // It may or may not be better, the dppd instruction is not particularly fast.

    return _mm_cvtsd_f64( res );
}

OpenMP 的 SIMD 指令可以向量化索引操作吗？

Can OpenMP's SIMD directive vectorize indexing operations?

c++

simd

openmp

eigen