如何在C/C++中为多线程FFTW启用OpenMP？代码是否线程安全？

在C/C++中实现FFTW时遭遇性能瓶颈，已按官方要求实现「仅创建一次FFTW计划、多次执行」的优化。现在尝试启用FFTW内置OpenMP并行功能，已用std::lock_guard<std::mutex>对计划创建等非线程安全操作加锁，但不确定fftw_plan_with_nthreads的调用时机，当前代码似乎未实现并行化。

代码示例：

#include"omp.h"  
#include <mutex>
#include <thread>
#include "fftw3.h"

static const int nx = 128;
static const int ny = 128;
static const int ncomp = 2;
static const int nyk = ny/2 + 1;

class MyFftwClass {
public:
    MyFftwClass(void);
    ~MyFftwClass(void);
    void execute(double rArr[], double cArr[][ncomp]);

private:
    static fftw_plan s_plan; // <-- shared by all instances
    double *m_buffer_in;
    fftw_complex *m_buffer_out;
};


class MyFftwClass1 {
public:
    MyFftwClass1(void);
    ~MyFftwClass1(void);
    void execute1(double cArr[][ncomp],double rArr[]);

private:
    static fftw_plan s_plan1; // <-- shared by all instances
    fftw_complex *m_buffer_in1;
    double *m_buffer_out1;
};

int main(){ 
    
fftw_init_threads(); //before calling any FFTW routines, you should call the function
//This function, which should only be called once (probably in your main() function), performs any one-time initialization required to use threads on your system.
int nThreads =1;//4;
omp_set_num_threads(nThreads);

double *Function;
Function= (double*) fftw_malloc(nx*ny*sizeof(double)); 
    for (int i = 0; i < nx; i++){
        for (int j = 0; j < ny; j++){
          Function[j + ny*i]  = //some initialization; 
            
        }
    }

fftw_complex *Functionk;
Functionk= (fftw_complex*) fftw_malloc(nx*nyk*sizeof(fftw_complex)); 
memset(Functionk, 42, nx*nyk* sizeof(fftw_complex));

MyFftwClass r2c1; //declare r2c1 of type MyFftwClass
r2c1.execute(Function,Functionk);

double *Out;
Out = (double*) fftw_malloc(nx*ny*sizeof(double)); 
memset(Out, 42, nx*ny* sizeof(double));

MyFftwClass1 c2r1; //declare r2c1 of type MyFftwClass
c2r1.execute1(Functionk,Out);
//fftw_free stuff
}
std::mutex g_my_fftw_mutex; // <-- must lock before using any fftw_*() function, except for fftw_execute*()
MyFftwClass1::MyFftwClass1(void)
{
    // serialize the initialization!
    std::lock_guard<std::mutex> lock(g_my_fftw_mutex);

    // allocate separate buffers for each instance
    m_buffer_in1 = fftw_alloc_complex(nx * nyk);
    m_buffer_out1 = fftw_alloc_real(nx * ny); 
    if (!(m_buffer_in1 && m_buffer_out1))
    {
        throw new std::runtime_error("Failed to allocate memory!");
    }

    // create plan *once* for all instances/threads
    if (!s_plan1)
    {
        s_plan1 = fftw_plan_dft_c2r_2d(nx, ny, m_buffer_in1, m_buffer_out1, FFTW_PATIENT);
        if (!s_plan1)
        {
            throw new std::runtime_error("Failed to create plan!");
        }
    }

}
MyFftwClass1::~MyFftwClass1(void)
{
    std::lock_guard<std::mutex> lock(g_my_fftw_mutex);
    fftw_destroy_plan(s_plan1);
    fftw_free(m_buffer_in1);
    fftw_free(m_buffer_out1);
}

void MyFftwClass1::execute1(double cArr[][ncomp],double rArr[]) //void MyFftwClass::execute(double *const in, fftw_complex *const out)
{
    // No serialization is needed here
    memcpy(m_buffer_in1,cArr, sizeof(fftw_complex) * nx*(nyk));
    fftw_execute_dft_c2r(s_plan1, m_buffer_in1, m_buffer_out1); //instead of fftw_excute(plan)
    memcpy(rArr, m_buffer_out1,sizeof(double) * nx*ny);
    //(rArr, 1.0 / (nx*ny), rArr); //renormalize 
}

fftw_plan MyFftwClass1::s_plan1 = NULL;

std::mutex g_my_fftw_mutex; // <-- must lock before using any fftw_*() function, except for fftw_execute*()
MyFftwClass::MyFftwClass(void)
//int r2cfft_initialize(r2cfft_t *const ctx, const std::size_t nx, const std::size_t ny)
{
    // serialize the initialization!
    std::lock_guard<std::mutex> lock(g_my_fftw_mutex);

    //allocating buffers first before plan creating gets rid off seg fault error

    // allocate separate buffers for each instance
    m_buffer_in = fftw_alloc_real(nx * ny); 
    m_buffer_out = fftw_alloc_complex(nx * nyk);
    if (!(m_buffer_in && m_buffer_out))
    {
        throw new std::runtime_error("Failed to allocate memory!");
    }

    // create plan *once* for all instances/threads
   if (!s_plan)
    {
        s_plan = fftw_plan_dft_r2c_2d(nx, ny, m_buffer_in, m_buffer_out, FFTW_PATIENT);
        if (!s_plan)
        {
            throw new std::runtime_error("Failed to create plan!");
        }
    }

    
}

MyFftwClass::~MyFftwClass(void)
{
    std::lock_guard<std::mutex> lock(g_my_fftw_mutex);
    fftw_free(m_buffer_in);
    fftw_free(m_buffer_out);
    fftw_destroy_plan(s_plan);
}

void MyFftwClass::execute(double rArr[], double cArr[][ncomp]) //void MyFftwClass::execute(double *const in, fftw_complex *const out)
{
    // No serialization is needed here
    memcpy(m_buffer_in, rArr,  sizeof(double) * nx*ny);
    fftw_execute_dft_r2c(s_plan, m_buffer_in, m_buffer_out); //instead of fftw_excute(plan)
    memcpy(cArr, m_buffer_out, sizeof(fftw_complex) * nx*(nyk));
}

fftw_plan MyFftwClass::s_plan = NULL;

请问这段代码是否线程安全？是否已实现并行化？

一、线程安全分析 #

当前代码存在线程安全隐患，主要问题如下：

1. 全局互斥量重复定义
   代码中两次定义了std::mutex g_my_fftw_mutex，会导致链接错误，全局变量不能重复定义。需将其改为单一全局变量，放在所有类定义之前。

2. 静态计划重复销毁
   两个类的静态计划s_plan/s_plan1会被每个实例的析构函数调用fftw_destroy_plan销毁。第一个实例析构时就会销毁计划，后续实例析构时会重复销毁已释放的计划，触发未定义行为。
   解决方法：为静态计划添加引用计数，仅当引用计数归0时才调用fftw_destroy_plan。

3. 正确的安全部分

   • 计划创建过程通过std::lock_guard加锁，避免多线程同时创建计划的冲突，符合FFTW线程安全规范。
   • 每个实例拥有独立的输入输出缓冲区，fftw_execute_dft_r2c/fftw_execute_dft_c2r是线程安全的（FFTW保证不同线程用不同缓冲区执行计划不会冲突），因此execute方法无需加锁，这部分逻辑正确。

二、并行化实现分析 #

当前代码未实现FFTW并行化，原因及修复方案如下：

1. 未调用fftw_plan_with_nthreads
   FFTW的并行计划需要在创建计划前设置线程数，正确调用顺序为：

   int main() {
       fftw_init_threads();
       int nThreads = 4; // 设置并行线程数
       fftw_plan_with_nthreads(nThreads); // 必须在创建计划前调用
       // ... 后续创建计划、执行FFT
   }
   

   该函数告诉FFTW在创建计划时生成并行化执行代码，缺少这一步，即使编译时启用OpenMP，计划也会是单线程的。

2. 错误使用omp_set_num_threads
omp_set_num_threads是OpenMP的API，对FFTW内置线程池无效。FFTW的并行线程数由fftw_plan_with_nthreads单独控制，无需调用omp_set_num_threads。

3. 无多线程测试场景
   main函数中仅创建单个实例执行FFT，即使计划是并行的，单线程执行也无法体现并行效果。需在多线程环境下（比如用std::thread创建多个线程同时调用execute）测试并行性能。

4. 编译需启用OpenMP支持
   除代码修改外，编译时要添加OpenMP编译选项（如GCC的-fopenmp），并链接FFTW并行版本库（如fftw3_omp）。

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