for loop parallelization sometimes works properly, not all t

Use this forum to discuss the book: Using OpenMP - Portable Shared Memory Parallel Programming, by Barbara Chapman, Gabriele Jost and Ruud van der Pas Read the http://openmp.org/forum/viewtopic.php?f=8&t=465 for book info and to download the examples. Post your feedback about the book and examples to this forum
Forum rules
The OpenMP Forums are now closed to new posts. Please visit Stack Overflow if you are in need of help: https://stackoverflow.com/questions/tagged/openmp
Locked
Arash_a

for loop parallelization sometimes works properly, not all t

Post by Arash_a »

Dear all,

Sorry, I just started learning OpenMP and there are lot of things to learn, so I am slightly confused.
I am analyzing my molecular dynamics simulation and in one part of the code I am trying to find the closest distance between water molecules (or ions) and protein. This is very time consuming part because I have around 500000 atoms and around 25000 frames. By single CPU it takes 1 week.
I changed this part of the code to parallel by OpenMP and it is really fast but with a little bug; 90% of the results (distances) are correct and 10% are wrong, compare to the single CPU code.
This is the part of my code which calculates the closest distance:

Code: Select all

 ...
for (i=0; i< number of frames(25000) )
...
// XP,YP,ZP protein coordinates; malloc allocation in the code
// XN,YN,ZN water molecules coordinates; malloc allocation
// LX,LY,LZ the dimension of simulation box, malloc allocation
// dimN defined as a temporary closest distance, malloc allocation
// NindN number of water molecules
…
...
l=0;kk=0;
#pragma omp parallel for shared(XN,YN,ZN,XP,YP,ZP,LX,LY,LZ,dimN,distN,xmin,ymin,zmin,i) private(kk,l)
      for (l=0; l < NindN; l++){

        if (XN[l]!=0.0 || YN[l]!=0.0 || ZN[l]!=0.0){

// this part relocates every thing inside a box with dimension LX*LY*LZ. xmin, ymin and zmin are the boundaries of the box.
        if (XN[l] < xmin) XN[l] += ceil((xmin - XN[l])/LX[i-1]) * LX[i-1];
        if (XN[l] > xmax) XN[l] -= ceil((XN[l] - xmax)/LX[i-1]) * LX[i-1];
        if (YN[l] < ymin) YN[l] += ceil((ymin - YN[l])/LY[i-1]) * LY[i-1];
        if (YN[l] > ymax) YN[l] -= ceil((YN[l] - ymax)/LY[i-1]) * LY[i-1];
        if (ZN[l] < zmin) ZN[l] += ceil((zmin - ZN[l])/LZ[i-1]) * LZ[i-1];
        if (ZN[l] > zmax) ZN[l] -= ceil((ZN[l] - zmax)/LZ[i-1]) * LZ[i-1];}

           for (kk=0; kk<NP; kk++){

         if (( XN[l]!=0.0 || YN[l]!=0.0 || ZN[l]!=0.0)   && ( XP[kk]!=0. || YP[kk]!=0. || ZP[kk]!=0. )  ){

         distN[l] = sqrt((XN[l]-XP[kk])*(XN[l]-XP[kk]) + (YN[l]-YP[kk])*(YN[l]-YP[kk]) + (ZN[l]-ZP[kk])*(ZN[l]-ZP[kk]) );
         if (distN[l] < dimN[l] ) {dimN[l] = distN[l];}

                  }
              }
            distN[l] = dimN[l];
           } 
#pragma omp barrier
...
Could you please tell me what is wrong with my code after parallelization? Why only in some cases it gives wrong answer and not all cases? I highly appreciate your comments and recommendation. I am using gcc on linux.
Thank you very much,

Cheers,
Arash

KentMilfeld

Re: for loop parallelization sometimes works properly, not a

Post by KentMilfeld »

Hi Arash,
I'm not sure why this fell through the cracks. Sorry about that. Did you resolve your problem. The updates to X/Y/ZN look straight forward to me, and the worksharing should not cause any problem. Are LX/Y/Z[i-1] defined for i=0? (At the end of your "for(l=0...) loop there is an implied barrier, so the "#pragma omp barrier" is redundant.)
Kent

gutha.raghu

Re: for loop parallelization sometimes works properly, not a

Post by gutha.raghu »

try this:

l=0;kk=0;
#pragma omp parallel for shared(XN,YN,ZN,XP,YP,ZP,LX,LY,LZ,dimN,distN,xmin,ymin,zmin,i) private(kk,l) schedule(static,100)
for (l=0; l < NindN; l++){

if (XN[l]!=0.0 || YN[l]!=0.0 || ZN[l]!=0.0){

// this part relocates every thing inside a box with dimension LX*LY*LZ. xmin, ymin and zmin are the boundaries of the box.
if (XN[l] < xmin) XN[l] += ceil((xmin - XN[l])/LX[i-1]) * LX[i-1];
if (XN[l] > xmax) XN[l] -= ceil((XN[l] - xmax)/LX[i-1]) * LX[i-1];
if (YN[l] < ymin) YN[l] += ceil((ymin - YN[l])/LY[i-1]) * LY[i-1];
if (YN[l] > ymax) YN[l] -= ceil((YN[l] - ymax)/LY[i-1]) * LY[i-1];
if (ZN[l] < zmin) ZN[l] += ceil((zmin - ZN[l])/LZ[i-1]) * LZ[i-1];
if (ZN[l] > zmax) ZN[l] -= ceil((ZN[l] - zmax)/LZ[i-1]) * LZ[i-1];}

for (kk=0; kk<NP; kk++){

if (( XN[l]!=0.0 || YN[l]!=0.0 || ZN[l]!=0.0) && ( XP[kk]!=0. || YP[kk]!=0. || ZP[kk]!=0. ) ){

distN[l] = sqrt((XN[l]-XP[kk])*(XN[l]-XP[kk]) + (YN[l]-YP[kk])*(YN[l]-YP[kk]) + (ZN[l]-ZP[kk])*(ZN[l]-ZP[kk]) );
if (distN[l] < dimN[l] ) {dimN[l] = distN[l];}

}
}
distN[l] = dimN[l];
}

...

Locked