使用 OpenMP 和 FFTW 时的不同结果答案

【问题标题】：Different results when using OpenMP and FFTW使用 OpenMP 和 FFTW 时的不同结果
【发布时间】：2016-10-03 15:41:17
【问题描述】：

我正在尝试并行化以下循环：

#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]) ;
        }
    }

我确实得到了一点加速，但无论我将哪些变量设置为 private 和 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

【问题讨论】：

为什么i+1 用于in 和out 的索引？
in 和 out 的大小为 2*Nx + 2，其余条目为 0。这样可以确保 in[0]=0
您是否构建/安装了 FFTW 的线程安全版本？请参阅 this answer，其中涵盖了有关将 FFTW 与 OpenMP 结合使用的此主题和其他相关主题。
我相信，因为我刚刚将函数 fftw_make_planner_thread_safe() 添加到该脚本中并且它没有给出编译错误。但它也没有给出正确的结果
什么是“p”？似乎可疑的是它没有在循环内修改，但是每次迭代都会调用一次 fftw_execute ，所以除非发生了什么神奇的事情，否则它肯定看起来好像你每次都在做同样的 FFT ......（他说不知道FFTW :-))。

标签： openmp fftw

【解决方案1】：

您在所有线程中使用相同的计划 p 调用 FFTW。由于计划包括输入和输出缓冲区的位置（提供给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);
}

【讨论】：

感谢您的回答！这导致“分段错误：11”
in 和 out 可能是指针而不是静态数组。在这种情况下，私有副本并不是真正私有的。在计划构造函数之前分配每个线程中的私有缓冲区，并在计划析构函数之后释放它们。
这导致另一个错误，我将发布一个编辑以显示我尝试过的内容
检查修改后的答案。如果它仍然不起作用，那么 fftw_execute 不是线程安全的，这是一个问题。另外，确保N 大于或等于Nx+1。
看来我遇到了更多麻烦，因为这也不起作用。我将尝试缩小范围并提出另一个问题。

【解决方案2】：

根本原因应该是这个patch没有向后移植到fftw-3.3.5版本，我认为你应该自己合并补丁。也可以参考讨论here。

【讨论】：