OpenMP并行循环内无法预创建FFTW计划的解决方案咨询

针对你在OpenMP并行for循环中多次执行FFTW变换的场景，这里提供几个可行的解决方案：

方案1：为每个线程创建私有FFTW计划 #

利用OpenMP的线程私有存储特性，让每个线程拥有独立的FFTW计划，避免共享数据冲突，同时复用计划提升性能。

实现步骤 #

• 用#pragma omp threadprivate声明计划变量，确保每个线程持有独立的计划实例
• 每个线程在首次执行迭代时创建对应大小的计划（因数据量固定为2的倍数，FFTW对此有优化）
• 每次迭代使用线程私有计划处理当前迭代的私有数据

代码示例 #

#include <fftw3.h>
#include <omp.h>

// 声明线程私有的FFTW计划
fftw_plan r2c_plan1, r2c_plan2, r2c_plan3, c2r_plan1, c2r_plan2;
#pragma omp threadprivate(r2c_plan1, r2c_plan2, r2c_plan3, c2r_plan1, c2r_plan2)

int main() {
    const int n = 1024; // 固定为2的倍数的数据量
    const int iterations = 1000;

    #pragma omp parallel for
    for (int i = 0; i < iterations; i++) {
        // 线程首次执行时初始化计划（需确保数据与计划生命周期匹配）
        if (!r2c_plan1) {
            double *in1 = fftw_malloc(sizeof(double) * n);
            fftw_complex *out1 = fftw_malloc(sizeof(fftw_complex) * (n/2 + 1));
            r2c_plan1 = fftw_plan_dft_r2c_1d(n, in1, out1, FFTW_ESTIMATE);
            // 同理初始化其他r2c和c2r计划
        }

        // 执行三次R2C变换
        fftw_execute_dft_r2c(r2c_plan1, private_in1, private_out1);
        fftw_execute_dft_r2c(r2c_plan2, private_in2, private_out2);
        fftw_execute_dft_r2c(r2c_plan3, private_in3, private_out3);

        // 处理变换后的复数数据（线程私有数据，无冲突）
        process_complex_data(private_out1, private_out2, private_out3, n);

        // 执行两次C2R变换
        fftw_execute_dft_c2r(c2r_plan1, private_complex_in1, private_real_out1);
        fftw_execute_dft_c2r(c2r_plan2, private_complex_in2, private_real_out2);
    }

    // 并行销毁线程私有计划并释放内存
    #pragma omp parallel
    {
        fftw_destroy_plan(r2c_plan1);
        fftw_destroy_plan(r2c_plan2);
        fftw_destroy_plan(r2c_plan3);
        fftw_destroy_plan(c2r_plan1);
        fftw_destroy_plan(c2r_plan2);
        // 释放对应的数据内存
        fftw_free(private_in1);
        // ... 其他数据释放操作
    }

    return 0;
}

方案2：每次迭代动态创建临时计划 #

如果线程私有计划的管理过于复杂，可以选择在每次迭代中创建临时计划，执行后立即销毁。由于数据量是2的倍数，FFTW_ESTIMATE标志创建计划的速度极快，1000次迭代的开销完全可接受。

代码示例 #

#include <fftw3.h>
#include <omp.h>

int main() {
    const int n = 1024;
    const int iterations = 1000;

    #pragma omp parallel for
    for (int i = 0; i < iterations; i++) {
        // 分配当前迭代的私有数据
        double *in1 = fftw_malloc(sizeof(double) * n);
        fftw_complex *out1 = fftw_malloc(sizeof(fftw_complex) * (n/2 + 1));
        // ... 分配其他变换所需数据

        // 1. 创建并执行第一次R2C变换
        fftw_plan r2c_p1 = fftw_plan_dft_r2c_1d(n, in1, out1, FFTW_ESTIMATE);
        fftw_execute(r2c_p1);
        fftw_destroy_plan(r2c_p1);

        // 2. 第二次R2C变换
        fftw_plan r2c_p2 = fftw_plan_dft_r2c_1d(n, in2, out2, FFTW_ESTIMATE);
        fftw_execute(r2c_p2);
        fftw_destroy_plan(r2c_p2);

        // 3. 第三次R2C变换
        fftw_plan r2c_p3 = fftw_plan_dft_r2c_1d(n, in3, out3, FFTW_ESTIMATE);
        fftw_execute(r2c_p3);
        fftw_destroy_plan(r2c_p3);

        // 4. 处理变换后的数据
        process_complex_data(out1, out2, out3, n);

        // 5. 第一次C2R变换
        fftw_plan c2r_p1 = fftw_plan_dft_c2r_1d(n, complex_in1, real_out1, FFTW_ESTIMATE);
        fftw_execute(c2r_p1);
        fftw_destroy_plan(c2r_p1);

        // 6. 第二次C2R变换
        fftw_plan c2r_p2 = fftw_plan_dft_c2r_1d(n, complex_in2, real_out2, FFTW_ESTIMATE);
        fftw_execute(c2r_p2);
        fftw_destroy_plan(c2r_p2);

        // 释放当前迭代的数据内存
        fftw_free(in1);
        fftw_free(out1);
        // ... 其他数据释放操作
    }

    return 0;
}

关键注意事项 #

• 所有变换数据必须是线程私有或每个迭代独立分配，绝对避免共享内存导致的竞争
• 优先选择方案1，复用线程私有计划能最大化FFTW的性能优势
• 若数据量在迭代中有变化（但你场景中固定为2的倍数），可在计划创建时动态传入当前数据大小

内容的提问来源于stack exchange，提问作者Patrick