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?

Answer 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).