Различных результатов при использовании OpenMP и FFTW

Question

Я пытаюсь распараллеливание следующего цикла:

#pragma omp parallel for private(j,i,mxy) firstprivate(in,out,p) 
    for(int j = 0; j < Ny; j++) { 
     //  #pragma omp parallel for private(i,mxy) firstprivate(in,my,j) 
     for(int i = 0; i < Nx; i++){ 
      mxy = i + j*Nx; 
      in[i+1] = b_2D[mxy] + I*0.0 ; 
     } 
     fftw_execute(p); 
     for(int i = 0; i < Nx; i++){ 
      mxy = i + j*Nx; 
      b_2D[mxy] = cimag(out[i+1]) ; 
     } 
    } 

Я получаю небольшую скорость вверх, но я получаю различный результат независимо от того, какого переменного я устанавливаю на частный и firstprivate. Я считаю, что это правильно, как я это сделал, но почему у меня получается другой результат, чем когда я запускаю это в серии?

Я попробовал следующее:

fftw_make_planner_thread_safe(); 
fftw_complex *in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); 
fftw_complex *out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); 

#pragma omp parallel private(j,i,mxy) firstprivate(in,out) 
{ 
    fftw_plan p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE); 
    for(j = 0; j < N; j++) 
     in[j] = 0.0; 

#pragma omp for 
    for(j = 0; j < Ny; j++) { 
     for(i = 0; i < Nx; i++) 
      in[i+1] = b_2D[i + j*Nx] + I*0.0; 
     fftw_execute(p); 
     for(i = 0; i < Nx; i++) 
      b_2D[i + j*Nx] = cimag(out[i+1]) ; 
    } 

    fftw_destroy_plan(p); 
} 
fftw_free(in); 
fftw_free(out); 

Это даст мне ошибку: "неисправность Сегментация: 11":

fftw_make_planner_thread_safe(); 
#pragma omp parallel private(j,i,mxy) 
{ 
    fftw_complex *in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); 
    fftw_complex *out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); 
    fftw_plan p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE); 
    for(j = 0; j < N; j++) 
     in[j] = 0.0; 

#pragma omp for 
    for(j = 0; j < Ny; j++) { 
     for(i = 0; i < Nx; i++) 
      in[i+1] = b_2D[i + j*Nx] + I*0.0; 
     fftw_execute(p); 
     for(i = 0; i < Nx; i++) 
      b_2D[i + j*Nx] = cimag(out[i+1]) ; 
    } 

    fftw_destroy_plan(p); 
    fftw_free(in); 
    fftw_free(out); 
} 

Я получаю эту ошибку снова

Если я это запустить : «Ошибка сегментации: 11» , но я снова забегаю, и он говорит:

solver(9674,0x7fff74e22000) malloc: *** error for object 0x7f8d70f00410: double free 
*** set a breakpoint in malloc_error_break to debug 
Abort trap: 6 

Answer 1

Вы звоните в FFTW с тем же планом p во все темы. Так как план включает расположение входных и выходных буферов (те, которые поставляются в конструктор плана fftw_plan_dft_whatever), все одновременные вызовы fftw_execute будут использовать те же буферы, а не частные копии. Решение состоит в том, чтобы построить отдельный план для каждого потока:

#pragma omp parallel private(j,i,mxy) firstprivate(in,out) 
{ 
    // The following OpenMP construct enforces thread-safety 
    // Remove if the plan constructor is thread-safe 
    #pragma omp critical (plan_ops) 
    fftw_plan my_p = fftw_plan_dft_whatever(..., in, out, ...); 
    // my_p now refers the private in and out arrays 

    #pragma omp for 
    for(int j = 0; j < Ny; j++) { 
     for(int i = 0; i < Nx; i++){ 
      mxy = i + j*Nx; 
      in[i+1] = b_2D[mxy] + I*0.0 ; 
     } 
     fftw_execute(my_p); 
     for(int i = 0; i < Nx; i++){ 
      mxy = i + j*Nx; 
      b_2D[mxy] = cimag(out[i+1]) ; 
     } 
    } 

    // See comment above for the constructor operation 
    #pragma omp critical (plan_ops) 
    fftw_destroy_plan(my_p); 
} 

Answer 2

первопричиной должно быть этим patch не портированные к FFTW-3.3.5 версии, и я думаю, вы должны объединить патч самостоятельно. Вы также можете обратиться к обсуждению here.