像这样使用 pragma omp simd 是否正确？

Question


#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
#define pow(x) ((x) * (x))
#define NUM_THREADS 8
#define wmax 1000
#define Nv 2
#define N 5

int b=0;

float Points[N][Nv]={ {0,1}, {3,4}, {1,2}, {5,1} ,{8,9}};

float length[wmax+1]={0};

float EuclDist(float* Ne, float* Pe) {
    int i;
    float s = 0;
 
    for (i = 0; i < Nv; i++) {
        s += pow(Ne[i] - Pe[i]); 
    }
    return s; 
}

void DistanceFinder(float* a[]){
    
    int i;
    #pragma omp simd
    for (i=1;i<N+1;i++){  
        length[b] += EuclDist(&a[i],&a[i-1]);
    }
    //printf(" %f\n", length[b]);
}



void NewRoute(){
//some irrelevant things

DistanceFinder(Points);
}

int main(){
    
omp_set_num_threads(NUM_THREADS);


do{
    b+=1;
    NewRoute();
    } while (b<wmax);
}

试图并行化这个循环并尝试不同的事情，尝试了这个。

似乎是最快的，但是这样使用 SIMD 是否正确？因为我使用的是之前的迭代（i 和 i - 1）。我看到的结果奇怪地正确与否。

Answer 1

Seems to be the fastest, however is it correct to use SIMD like that?

首先，有一个竞争条件需要修复，即在数组更新期间length[b]。此外，您正在访问数组 a 之外的内存；（从 1 迭代到 N + 1），你正在传递 &a[i]。您可以使用 OpenMP reduction 子句修复 竞争条件 ：

void DistanceFinder(float* a[]){
    
    int i;
    float sum = 0;
    float tmp;
    #pragma omp simd private(tmp) reduction(+:sum)
    for (i=1;i<N;i++){  
        tmp = EuclDist(a[i], a[i-1]);
        sum += tmp;
    }
    length[b] += sum;
}

另外需要提供EuclDist的版本如下：

#pragma omp declare simd uniform(Ne, Pe)
float EuclDist(float* Ne, float* Pe) {
    int i;
    float s = 0;
    for (i = 0; i < Nv; i++)
        s += pow(Ne[i] - Pe[i]); 
    return s; 
}

Because I'm using a previous iteration (i and i - 1).

在你的情况下，没关系，因为数组 a 刚刚被读取。

The results I see though are correct weirdly or not.

很可能没有进行矢量化。无论如何，由于前面提到的竞争条件。

，它仍然是未定义的行为

您可以简化代码，以增加矢量化实际发生的可能性，例如：

void DistanceFinder(float* a[]){
    int i;
    float sum = 0;
    float tmp;
    #pragma omp simd private(tmp) reduction(+:sum)
    for (i=1;i<N;i++){  
        tmp = pow(a[i][0] - a[i-1][0]) + pow(a[i][1] - a[i-1][1]) 
        sum += tmp;
    }
    length[b] += sum;
}

为了提高代码的性能，您可以做的进一步更改是分配矩阵（作为函数 DistanceFinder 的参数传递），其方式是在迭代其行时（即 a[i]）您将迭代连续的内存地址。

例如，您可以传递两个数组 a1 和 a2 来表示矩阵的第一列和第二列 a:

  void DistanceFinder(float a1[], float a2[]){
        int i;
        float sum = 0;
        float tmp;
        #pragma omp simd private(tmp) reduction(+:sum)
        for (i=1;i<N;i++){  
            tmp = pow(a1[i] - a1[i-1]) + pow(a2[i][1] - a2[i-1][1]) 
            sum += tmp;
        }
        length[b] += sum;
    }

像这样使用 pragma omp simd 是否正确？

Is using pragma omp simd like this correct?

c

parallel-processing

simd

openmp