Векторизация на непоследовательных элементах через набор индексов

Question

Стандартный шаблон для векторизации выглядит так:

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

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

, где все индексы доступны последовательно.

Однако у меня есть ситуация, когда не все N элементы arr нуждаются в обновлении. Шаблон у меня выглядит следующим образом:

#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.

В частности, какие из приведенных выше рекомендаций гарантируют безопасную векторизацию?

Answer 1

На основании дополнительных комментариев предполагается, что функция, которую вы хотите оптимизировать, выглядит следующим образом:

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);
    }
}

Поскольку ваша цель - Pentium, доступны только инструкции SSE2 SIMD. Вы можете попробовать следующую оптимизированную функцию (также доступна здесь ):

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

Вы могли бы также выполнить некоторое развертывание цикла (это, вероятно, сделал бы векторизатор компилятора). Но в этом случае тело цикла довольно существенно и, вероятно, связано с вводом / выводом памяти, поэтому я не думаю, что это сильно поможет.

Если значение этих функций неясно, вы можете прочитать this .

Всестороннее обсуждение векторизации здесь . Короче говоря, современные процессоры предлагают векторные регистры в дополнение к обычным. В зависимости от архитектуры машины у нас есть 128-битные регистры (инструкции SSE), 256-битные регистры (инструкции AVX), 512-битные регистры (инструкции AVX512). Векторизация означает использование этих регистров в их полном объеме. Например, если у нас есть вектор double {x0, x1, x2, …. xn}, и мы хотим добавить его к константе k, получая {x0+k, x1+k, x2+k, …. xn+k}, в скалярном режиме операция разлагается на чтение x(i), добавление k, сохранить результат. В векторном режиме с SSE это выглядело бы как чтение x[i] и x[i+1] одновременно, добавление k одновременно, сохранение 2 результатов одновременно. Если элементы не являются смежными в памяти, мы не можем загружать x [i] и x [i + 1] одновременно, а также не можем сохранять результаты одновременно, но, по крайней мере, мы можем выполнить два сложения одновременно.