关于OpenMP和-Ofast的结合

Question

我在 for 循环中实现了 OpenMP 并行化，其中我有一个 sum，这是导致我的代码变慢的主要原因。当我这样做时，我发现最终结果与我为非并行化代码（用 C 编写）获得的结果不同。所以首先，人们可能会想“好吧，我只是没有很好地实现并行化”，但奇怪的是，当我使用 -Ofast 优化运行并行化代码时，结果突然是正确的。

那就是：

• -O0 正确
• -Ofast 正确
• OMP -O0 错误
• OMP -O1 错误
• OMP -O2 错误
• OMP -O3 错误
• OMP -Ofast 正确！

-Ofast 可以做什么来解决仅在我实施 openmp 时出现的错误？ 关于我可以检查或测试什么的任何建议？ 谢谢！

编辑 在这里，我包含了仍然重现问题的最小版本的代码。

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

#define LENGTH 100
#define R 50.0
#define URD 1.0/sqrt(2.0)
#define PI (4.0*atan(1.0)) //pi

const gsl_rng_type * Type;
gsl_rng * item;

double CalcDeltaEnergy(double **M,int sx,int sy){
    double DEnergy,r,zz;
    int k,j;
    double rrx,rry;
    int rx,ry;
    double Energy, Cpm, Cmm, Cmp, Cpp;
    DEnergy = 0;

    //OpenMP parallelization:
    #pragma omp parallel for reduction (+:DEnergy)
    for (int index = 0; index < LENGTH*LENGTH; index++){
        k = index % LENGTH;
        j = index / LENGTH;

    zz = 0.5*(1.0 - pow(-1.0, k + j + sx + sy));
    for (rx = -1; rx <= 1; rx++){
        for (ry = -1; ry <= 1; ry++){
            rrx = (sx - k - rx*LENGTH)*URD;
            rry = (sy - j - ry*LENGTH)*URD;

            r = sqrt(rrx*rrx + rry*rry + zz);
            if(r != 0 && r <= R){
                Cpm = sqrt((rrx+0.5*(0.702*cos(M[k][j])-0.702*cos(M[sx][sy])))*(rrx+0.5*(0.702*cos(M[k][j])-0.702*cos(M[sx][sy]))) + (rry+0.5*(0.702*sin(M[k][j])-0.702*sin(M[sx][sy])))*(rry+0.5*(0.702*sin(M[k][j])-0.702*sin(M[sx][sy]))) + zz);
                Cmm = sqrt((rrx-0.5*(0.702*cos(M[k][j])-0.702*cos(M[sx][sy])))*(rrx-0.5*(0.702*cos(M[k][j])-0.702*cos(M[sx][sy]))) + (rry-0.5*(0.702*sin(M[k][j])-0.702*sin(M[sx][sy])))*(rry-0.5*(0.702*sin(M[k][j])-0.702*sin(M[sx][sy]))) + zz);
                Cpp = sqrt((rrx+0.5*(0.702*cos(M[k][j])+0.702*cos(M[sx][sy])))*(rrx+0.5*(0.702*cos(M[k][j])+0.702*cos(M[sx][sy]))) + (rry+0.5*(0.702*sin(M[k][j])+0.702*sin(M[sx][sy])))*(rry+0.5*(0.702*sin(M[k][j])+0.702*sin(M[sx][sy]))) + zz);
                Cmp = sqrt((rrx-0.5*(0.702*cos(M[k][j])+0.702*cos(M[sx][sy])))*(rrx-0.5*(0.702*cos(M[k][j])+0.702*cos(M[sx][sy]))) + (rry-0.5*(0.702*sin(M[k][j])+0.702*sin(M[sx][sy])))*(rry-0.5*(0.702*sin(M[k][j])+0.702*sin(M[sx][sy]))) + zz);
                Cpm = 1.0/Cpm;
                Cmm = 1.0/Cmm;
                Cpp = 1.0/Cpp;
                Cmp = 1.0/Cmp;
                Energy = (Cpm + Cmm - Cpp - Cmp)/(0.702*0.702); // S=cte=1

                DEnergy -= 2.0*Energy;
            }
        }
    }
    }
return DEnergy;
}

void Initialize(double **M){
double random;
    for(int i=0;i<(LENGTH-1);i=i+2){
          for(int j=0;j<(LENGTH-1);j=j+2) {
              random=gsl_rng_uniform(item);
              if (random<0.5) M[i][j]=PI/4.0;
              else M[i][j]=5.0*PI/4.0;

              random=gsl_rng_uniform(item);
              if (random<0.5) M[i][j+1]=3.0*PI/4.0;
              else M[i][j+1]=7.0*PI/4.0;

              random=gsl_rng_uniform(item);
              if (random<0.5) M[i+1][j]=3.0*PI/4.0;
              else M[i+1][j]=7.0*PI/4.0;

              random=gsl_rng_uniform(item);
              if (random<0.5) M[i+1][j+1]=PI/4.0;
              else M[i+1][j+1]=5.0*PI/4.0;
          }
    }
}

int main(){
    //Choose and initiaze the random number generator
    gsl_rng_env_setup();
    Type = gsl_rng_default; //default=mt19937, ran2, lxs0
    item = gsl_rng_alloc (Type);

    double **S; //site matrix
    S = (double **) malloc(LENGTH*sizeof(double *));
    for (int i = 0; i < LENGTH; i++)
        S[i] = (double *) malloc(LENGTH*sizeof(double ));

    //Initialization
    Initialize(S);

    int l,m;
    for (int cl = 0; cl < LENGTH*LENGTH; cl++) {
        l = gsl_rng_uniform_int(item, LENGTH); // RNG[0, LENGTH-1]
        m = gsl_rng_uniform_int(item, LENGTH); // RNG[0, LENGTH-1]
        printf("%lf\n", CalcDeltaEnergy(S, l, m));
    }


    //Free memory
    for (int i = 0; i < LENGTH; i++)
        free(S[i]);
    free(S);
    return 0;
} 

我编译：

g++ [optimization] -lm test.c -o test.x -lgsl -lgslcblas -fopenmp

并运行：

GSL_RNG_SEED=123; ./test.x > test.dat

比较不同优化的输出可以看到我之前所说的。

Answer 1

免责声明：我几乎没有使用 OpenMP 的经验

这可能是您在使用 OpenMP 时遇到的一种竞争条件。

您需要将 OpenMP 循环内的所有这些变量声明为私有。一个内核可能会针对某个值 index 计算它们的值，该值会在使用另一个值 index 的内核上迅速重新计算为不同的值：变量如 k、j、rrx、 rry 等在计算节点之间共享。

而不是使用类似的编译指示

#pragma omp parallel for private(k,j,zz,rx,ry,rrx,rry,r,Cpm,Cmm,Cpp,Cmp,Energy) reduction (+:D\

（以下 Zulan 评论的学分：）您还可以在并行区域内声明变量，尽可能在本地。这使它们隐式私有化，不易出现初始化问题且更易于推理。

（您甚至可以考虑将所有内容都放在外部 for 循环中（超过 index）在一个函数中：与计算相比，函数调用开销是最小的。）

至于为什么-Ofast 与 OpenMP 一起确实会产生正确的输出。

我的猜测是：主要是运气。这是-Ofast 所做的（gcc 手册）：

无视严格的标准合规性。 -Ofast 启用所有 -O3 优化。它还支持并非对所有符合标准的程序都有效的优化。它打开 -ffast-math [...]

这是-ffast-math的部分：

除 -Ofast 之外的任何 -O 选项都不会打开此选项，因为它可能导致依赖于数学函数的 IEEE 或 ISO 规则/规范的精确实现的程序的错误输出。但是，对于不需要这些规范保证的程序，它可能会产生更快的代码。

因此，sqrt、cos 和 sin 可能会更快。我的猜测是，在这种情况下，外循环内变量的计算不会相互咬合，因为各个线程是如此之快，它们不会发生冲突。但这是一个非常随意的解释和猜测。