Making functions thread safe in C for OpenMP

Question

I have a function which determines if an array needs to be smoothed. I have multiple arrays that this operation needs to be preformed on. I'd like to use the sections construct in OpenMP to do this. Unfortunately, I'm running into segmentation faults as soon as the code grows in an appreciable size. This may be a memory limitation issue but I'd like your input. Here's the pseudo code for the call to the function:

#pragma omp parallel default(shared) private(i,j,k,memd,memi,mdpmi,mdpme,mipmn,mipmn,mdinv,meinv,miinv,mimime,memime,gchk,spat,dtinv,mtempx,mtempy,mtempz,mtempinv) \
                     firstprivate(memd,memi,mdpmi,mdpme,mipmn,mipmn,mdinv,meinv,miinv,mimime,memime,gchk,spat,dtinv)
{
.
.
.
  #pragma omp single
  printf("----- Check for smoothing -----\n");
  #pragma omp sections
  {
    //Grids are now at the same timestep so we smooth
    #pragma omp section
    {
      printf("-----   Smoothing of Rho by TID:%d\n",omp_get_thread_num());
      smooth(mhd->rho,mhd->rhosmo,grid,omp_get_thread_num());
    }
    #pragma omp section
    {
      printf("-----   Smoothing of Rhoi by TID:%d\n",omp_get_thread_num());
      smooth(mhd->rhoi,mhd->rhoismo,grid,omp_get_thread_num());
    }
    .
    .
    .
  } /*-- End of Sections --*/
  .
  .
  .
} /*-- End of Parallel Region --*/

Now the function looks like this

void smooth(double ***field,char ***tsmooth,GRID *grid,int tid)
{
  double mtempx[grid->nx][grid->ny][grid->nz];
  double mtempy[grid->nx][grid->ny][grid->nz];
  double mtempz[grid->nx][grid->ny][grid->nz];
  double mtempinv;
  double etol=1e-10;  //Oscillation amplitude tollerance
  double itol=1e-2;   //Inverse tollerance
  for(i=0;i<grid->nx;i++)
  {
    for(j=0;j<grid->ny;j++)
    {
      for(k=0;k<grid->nz;k++)
        {
          printf("-----  SMOOTH(TID:%2d) i=%3d j=%3d k=%3d\n",tid,i,j,k);
          mtempx[i][j][k]=0.0;
          mtempy[i][j][k]=0.0;
          mtempz[i][j][k]=0.0;
          tsmooth[i][j][k]=0;
        }
    }
  }
  for(i=1;i<grid->nx-1;i++)
  {
    for(j=1;j<grid->ny-1;j++)
    {
      for(k=1;k<grid->nz-1;k++)
      {
        mtempinv=1.;
        if (sqrt(field[i][j][k]*field[i][j][k]) > itol) mtempinv=1./field[i][j][k];
        mtempx[i][j][k]=(field[i+1][j][k]-field[i][j][k])*(field[i][j][k]-field[i-1][j][k]);
        mtempy[i][j][k]=(field[i][j+1][k]-field[i][j][k])*(field[i][j][k]-field[i][j-1][k]);
        mtempz[i][j][k]=(field[i][j][k+1]-field[i][j][k])*(field[i][j][k]-field[i][j][k-1]);
        if (sqrt(mtempx[i][j][k]*mtempx[i][j][k])*mtempinv*mtempinv <= etol) mtempx[i][j][k]=0.0;
        if (sqrt(mtempy[i][j][k]*mtempy[i][j][k])*mtempinv*mtempinv <= etol) mtempy[i][j][k]=0.0;
        if (sqrt(mtempz[i][j][k]*mtempz[i][j][k])*mtempinv*mtempinv <= etol) mtempz[i][j][k]=0.0;
      }
    }
  }
  for(i=1;i<grid->nx-1;i++)
  {
    for(j=1;j<grid->ny-1;j++)
    {
      for(k=1;k<grid->nz-1;k++)
      {
        if (   ((mtempx[i][j][k] < 0.0) && ((mtempx[i+1][j][k] < 0.0)||(mtempx[i-1][j][k] < 0.0))) 
                || ((mtempy[i][j][k] < 0.0) && ((mtempy[i][j+1][k] < 0.0)||(mtempy[i][j-1][k] < 0.0)))
                || ((mtempz[i][j][k] < 0.0) && ((mtempz[i][j][k+1] < 0.0)||(mtempx[i][j][k-1] < 0.0)))) tsmooth[i][j][k]=1;
      }
    }
  }

}

Now this worked fine in serial and for arrays 101x7x49 it seemed to work fine but when I went to an array of 389x7x739 I could get past the first printf statement in the function. Note I only added the printf statements to help understand what was going on. Shouldn't this be threadsafe?

Throwback1986 · Accepted Answer · 2010-10-26 18:12:56Z

1

The temporary matrices (mtempx, mtempy, etc.) require more space than your system stack provides. Dynamically allocate those buffers instead (or statically declare them to be a fixed size).

answered Oct 26, 2010 at 18:12

Throwback1986

5,9951 gold badge35 silver badges23 bronze badges

Sign up to request clarification or add additional context in comments.

3 Comments

Lazer Over a year ago

I believe you are correct. I'm going to try dynamically allocating them.

Throwback1986 Over a year ago

And, to answer the originally posed question: I don't see any obvious thread safety issues in the function.

CocoCrisp Over a year ago

If one stores the matrix in data segment of memory it solves the problem as data segment has more memory than stack frame. Are there any pitfalls to it except that the data frame becomes full as well?

Collectives™ on Stack Overflow

Making functions thread safe in C for OpenMP

1 Answer 1

3 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

3 Comments

Your Answer

Sign up or log in

Post as a guest

Related