2010-10-26 71 views
0

我有一個函數可以確定數組是否需要平滑。我有多個陣列需要執行此操作。我想使用OpenMP中的部分構造來完成此操作。不幸的是,只要代碼增長得很快,我就會遇到分段錯誤。這可能是一個內存限制問題,但我希望你的意見。 下面是調用該函數的僞代碼:在C中爲OpenMP使線程安全函數安全

#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 --*/ 

現在的功能看起來像這樣

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; 
     } 
    } 
    } 

}

現在這個串行和數組101x7x49它似乎工作得很好工作正常,但當我去了一個389x7x739的數組時,我可以通過函數中的第一個printf語句。注意我只添加了printf語句來幫助理解發生了什麼。這不應該是線程安全嗎?

回答

1

臨時矩陣(mtempx,mtempy等)需要比系統堆棧提供的更多空間。改爲動態分配這些緩衝區(或靜態聲明它們是固定大小)。

+0

我相信你是對的。我將嘗試動態分配它們。 – Lazer 2010-10-26 19:08:01

+0

而且,爲了回答最初提出的問題:我沒有在函數中看到任何明顯的線程安全問題。 – Throwback1986 2010-10-26 19:39:06