c - OpenMP for 循环没有给出正确的结果

Question

抱歉，我刚开始学习 OpenMP，所以有点困惑。我正在分析我的分子动力学模拟，在代码的一部分中，我试图找到水分子(或 ionic )和蛋白质之间的最近距离。这是非常耗时的部分，因为我有大约 500000 个原子和大约 25000 帧。单 CPU 需要 1 周时间(一组计算不仅是距离)。

我将这部分代码更改为通过 OpenMP 进行并行处理，速度非常快，但有一点错误；与单 CPU 代码相比，90% 的结果(距离)是正确的，10% 是错误的。这是我的代码中计算最近距离的部分:

    ...
    for (i=0; i< number of frames(25000) )
    ...
    // XP,YP,ZP protein coordinates; malloc allocation in the code
    // XI,YI,ZI Sodium molecule coordinates; malloc allocation
    // LX,LY,LZ the dimension of simulation box, malloc allocation
    // dimI defined as a temporary closest distance, filled with very large constant,
    // malloc allocation
    // NSOD number of Sodium molecules
    // rhos keeping the closest distance for each Sodium for each frame.
    …
    ...int l=0,kk=0;
#pragma omp parallel for shared(XI,YI,ZI,XP,YP,ZP,LX,LY,LZ,qq,dimI,distI,rhos,xmin,ymin,zmin,i) private(kk,l)
      for (l=0; l < NSOD; l++){
        // this part relocates every thing inside a box with dimension LX*LY*LZ. xmin, ymin and zmin are the boundaries of the box.
        if (XI[l]!=0.0 || YI[l]!=0.0 || ZI[l]!=0.0){
          if (XI[l] < xmin) XI[l] += ceil((xmin - XI[l])/LX[i-1]) * LX[i-1];
          if (XI[l] > xmax) XI[l] -= ceil((XI[l] - xmax)/LX[i-1]) * LX[i-1];
          if (YI[l] < ymin) YI[l] += ceil((ymin - YI[l])/LY[i-1]) * LY[i-1];
          if (YI[l] > ymax) YI[l] -= ceil((YI[l] - ymax)/LY[i-1]) * LY[i-1];
          if (ZI[l] < zmin) ZI[l] += ceil((zmin - ZI[l])/LZ[i-1]) * LZ[i-1];
          if (ZI[l] > zmax) ZI[l] -= ceil((ZI[l] - zmax)/LZ[i-1]) * LZ[i-1];
        }
        for (kk=0; kk<NP; kk++){

         if ( ( XP[kk]!=0. || YP[kk]!=0. || ZP[kk]!=0. )  ){
           distI[l] = sqrt((XI[l]-XP[kk])*(XI[l]-XP[kk]) + (YI[l]-YP[kk])*(YI[l]-YP[kk]) + (ZI[l]-ZP[kk])*(ZI[l]-ZP[kk]) );
           if (distI[l] < dimI[l] ) {
              dimI[l] = distI[l];
           }
         }
        }
        distI[l] = dimI[l];
        rhos[qq][l] = dimI[l];

} #pragma omp 屏障 ...

你能告诉我并行化后我的代码有什么问题吗？为什么只有在某些情况下它会给出错误的答案而不是在所有情况下？我非常感谢您的意见和建议。我在 Linux 上使用 gcc。非常感谢，

干杯，阿拉什

Answer 1

在处理 float 时，拥有

可能不是一个好主意

if (XI[l]!=0.0 || YI[l]!=0.0 || ZI[l]!=0.0){

相反，您应该与 epsilon(是一些非常小的数字)进行比较

if (fabs(XI[l]) > epsilon || ...

否则这可能会导致问题。