c - 关于 OpenMP 和 -Ofast 的组合

Question

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

那就是:

• -O0 正确
• -快速正确
• 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 循环内的所有这些变量声明为私有(private)。一个核心可能会根据 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 的评论:)您还可以尽可能在本地声明并行区域内的变量。这使它们隐式私有(private)并且不太容易出现初始化问题并且更容易推理。

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

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

我的猜测是:主要是运气。以下是 -Ofast 的作用(gcc 手册):

Disregard strict standards compliance. -Ofast enables all -O3 optimizations. It also enables optimizations that are not valid for all standard-compliant programs. It turns on -ffast-math [...]

这是关于 -ffast-math 的部分:

This option is not turned on by any -O option besides -Ofast since it can result in incorrect output for programs that depend on an exact implementation of IEEE or ISO rules/specifications for math functions. It may, however, yield faster code for programs that do not require the guarantees of these specifications.

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