Myujth

The standard template for vectorization seems to be thus:

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

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

where all of the indices are consecutively accessed.

However, I have a situation where not all Nelements of arr need updation. The template I have is as follows:

#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]);

In this case, Intel Advisor reports that this loop was not vectorized because loop iteration count cannot be computed before executing the loop. In my application, this loop is executed for different subsets of indices, and for each such subset, number_in_index_set and indexset would change correspondingly.

I have two questions:

(1)What does it mean for the problematic loop to even be vectorized? The array indices are not consecutive, so how would the compiler even go about vectorizing the loop?

(2)Assuming vectorization is possible, as Intel Advisor seems to suggest, how can the loop in question be safely vectorized? The recommendation from Intel Advisor are thus:

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.

Specifically, which of the above recommendations are guaranteed to ensure safe vectorization?

Based on the additional comments, assuming the function you want to optimize is the following:

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

Since your target is Pentium, only SSE2 SIMD instructions are available. You can try the following optimized function (also available 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);
 

Here the load and store operations are not vectorized, i.e. we are loading x and y one by one and we are storing into result one by one. All other operations are vectorized.

With gcc the assembler output of the body of the main loop is below.

 // 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

You could also do some loop unrolling (this is probably what a compiler vectorizer would do). But in this case the loop body is quite substantial, and probably memory I/O bound, so I do not think it would help much.

If the meaning of these functions is not clear, you can read this.

A comprehensive discussion on vectorization is here. In brief, modern cpu offer vectorial registers I addition to the regular ones. Depending on the machine architecture, we have 128 bit registers (SSE instructions), 256 bit registers (AVX instructions), 512 bit registers (AVX512 instructions). Vectorizing means using these register at their full capability. For instance if we have a vector of double x0, x1, x2, …. xn and we want to add it to a constant k, obtaining x0+k, x1+k, x2+k, …. xn+k, in scalar mode the operation is decomposed in read x(i), add k, store result. In vectorial mode, with SSE, it would look like read x[i] and x[i+1] simultaneously, add k simultaneously, store 2 results simultaneously.
If the elements are not contiguous in memory, we cannot load x[i] and x[i+1] simultaneously, nor we can store the results simultaneously, but at least we can perform the two additions simultaneously.