矢量化的标准模板似乎是这样的:

#define N 100
double arr[N];
double func(int i);

for(int i = 0; i <N; i++)
    arr[i] = func(i);

连续访问所有索引的位置。

但是，我遇到的情况不是 arr 的所有 N 元素都需要更新。我的模板如下:

#define N 100
double arr[N];
double func(int i);

int indexset[N];//this indexset has the indices of arr[] that get updated
int number_in_index_set;
//E.g., if I only need to update arr[4] and arr[10], number_in_index_set = 2
//indexset[0] = 4 and indexset[1] = 10

for(int i = 0; i <number_in_index_set; i++)
    arr[indexset[i]] = func(indexset[i]);

在这种情况下，Intel Advisor 报告此循环未矢量化，因为无法在执行循环之前计算循环迭代次数。在我的应用程序中，这个循环是针对不同的索引子集执行的，对于每个这样的子集，number_in_index_set 和 indexset[] 都会相应地改变。

我有两个问题:

(1)有问题的循环甚至被矢量化意味着什么？数组索引不是连续的，那么编译器将如何对循环进行矢量化？

(2) 假设矢量化是可能的，正如 Intel Advisor 似乎暗示的那样，如何才能安全地将所讨论的循环矢量化？因此，Intel Advisor 的建议是:

For Example 1, where the loop iteration count is not available before the loop 

executes: If the loop iteration count and iterations lower bound can be calculated 

for the whole loop:
    Move the calculation outside the loop using an additional variable.
    Rewrite the loop to avoid goto statements or other early exits from the loop.
    Identify the loop iterations lower bound using a constant.
For example, introduce the new limit variable:
void foo(float *A) {
  int i;
  int OuterCount = 90;
  int limit = bar(int(A[0]));
  while (OuterCount > 0) {
    for (i = 1; i < limit; i++) {
      A[i] = i + 4;
    }
    OuterCount--;
  }
}
For Example 2, where the compiler cannot determine if there is aliasing between all 

the pointers used inside the loop and loop boundaries: Assign the loop boundary value 

to a local variable. In most cases, this is enough for the compiler to determine 

aliasing may not occur.
You can use a directive to accomplish the same thing automatically.
Target  ICL/ICC/ICPC Directive
Source loop     #pragma simd or #pragma omp simd
Do not use global variables or indirect accesses as loop boundaries unless you also 

use one of the following:
    Directive to ignore vector dependencies
    Target  ICL/ICC/ICPC Directive
    Source loop     #pragma ivdep
    restrict
    keyword.

具体来说，上述哪些建议可以保证安全矢量化？

最佳答案

根据补充的意见，假设你要优化的功能如下:

void euclidean(double x0, double y0, const double *x, const double* y,
               const size_t *index, size_t index_size, double *result)
{
    for (size_t i = 0; i < index_size; ++i) {
        double dx = x0 - x[index[i]];
        double dy = y0 - y[index[i]];
        result[index[i]] = sqrt(dx*dx + dy * dy);
    }
}

由于您的目标是奔腾，因此只有 SSE2 SIMD 指令可用。您可以尝试以下优化功能(也可用 here ):

void euclidean_sse2(double x0, double y0, const double *x, const double* y,
               const size_t *index, size_t index_size, double *result)
{
    __m128d vx0 = _mm_set1_pd(x0);
    __m128d vy0 = _mm_set1_pd(y0);
    for (size_t i = 0; i + 1 < index_size; i += 2) { // process 2 at a time
        size_t i0 = index[i];
        size_t i1 = index[i+1];
        __m128d vx = _mm_set_pd(x[i1], x[i0]);  // load 2 scattered elements
        __m128d vy = _mm_set_pd(y[i1], y[i0]);  // load 2 scattered elements
        __m128d vdx = _mm_sub_pd(vx, vx0);
        __m128d vdy = _mm_sub_pd(vy, vy0);
        __m128d vr = _mm_sqrt_pd(_mm_add_pd(_mm_mul_pd(vdx, vdx), _mm_mul_pd(vdy, vdy)));
        _mm_store_sd(result + i0, vr);
        _mm_storeh_pd(result + i1, vr);
    }
    if (index_size & 1) {  // if index_size is odd, there is one more element to process
        size_t i0 = index[index_size-1];
        double dx = x0 - x[i0];
        double dy = y0 - y[i0];
        result[i0] = sqrt(dx*dx + dy * dy);
    }
}

这里的加载和存储操作不是向量化的，即我们正在一个一个地加载 x 和 y 并且我们正在存储到 result逐个。所有其他操作都是向量化的。

对于 gcc，主循环体的汇编器输出如下。

    // scalar load operations
    mov     r10, QWORD PTR [rax]
    mov     r9, QWORD PTR [rax+8]
    add     rax, 16
    movsd   xmm3, QWORD PTR [rdi+r10*8]
    movsd   xmm2, QWORD PTR [rsi+r10*8]
    movhpd  xmm3, QWORD PTR [rdi+r9*8]
    movhpd  xmm2, QWORD PTR [rsi+r9*8]
    // vectorial operations
    subpd   xmm3, xmm4
    subpd   xmm2, xmm5
    mulpd   xmm3, xmm3
    mulpd   xmm2, xmm2
    addpd   xmm2, xmm3
    sqrtpd  xmm2, xmm2
    // scalar store operations
    movlpd  QWORD PTR [r8+r10*8], xmm2
    movhpd  QWORD PTR [r8+r9*8], xmm2

您还可以进行一些循环展开(这可能是编译器向量化器会做的)。但在这种情况下，循环体相当大，并且可能受内存 I/O 限制，所以我认为它不会有太大帮助。

如果对这些函数的含义不清楚，可以阅读this .

关于矢量化的全面讨论是 here .简而言之，现代 cpu 除了常规寄存器外还提供 vector 寄存器。根据机器架构的不同，我们有 128 位寄存器(SSE 指令)、256 位寄存器(AVX 指令)、512 位寄存器(AVX512 指令)。矢量化意味着充分利用这些寄存器。例如，如果我们有一个 vector double {x0, x1, x2, …. xn}，我们想把它加到常量k上，得到{x0+k, x1+k, x2+k, …. xn+k}，在标量模式下，操作分解为读取x(i)，添加k，存储结果。在 vector 模式下，使用 SSE，看起来像是同时读取 x[i] 和 x[i+1]，同时添加 k，同时存储 2 个结果。 如果元素在内存中不连续，我们不能同时加载x[i]和x[i+1]，也不能同时存储结果，但至少我们可以同时进行这两个加法。

关于c++ - 通过索引集对非连续元素进行矢量化，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/53246116/