使用 std::vector 容器进行 FFT

问题描述

我想请教您关于在执行 FFT 中使用 std::vector 容器的建议。

在我的程序中，我在循环中使用某个函数 (get_PDF_CDF)。当我使用 std::vector 容器执行 FFT 时，在某个迭代（通常是第二次）期间，我收到了分段错误/内存损坏错误。

我知道 FFTW 最好与 fftw_malloc 和 fftw_free 提供的对齐配合使用，但我想增加一些趣味并在整个程序中始终如一地使用 std::vector。

我已经重写了函数以使用推荐的分配例程 - 它基本上将向量的内容复制到使用 fftw 例程分配和释放的浮点数组 - ~~并且它完美地工作~~ 但它仍然有问题。它还涉及对 gaussian_f 和 top_hat_f 过滤器的简单重写。此外，我想了解是否可以使用 std::vectors 以任何方式修改版本，使其不会产生分段错误。

以下是两个版本，截取其基本内容。如果这不足以进行错误跟踪，我深表歉意 - 在这种情况下，我将提供一个工作示例，尽管从代码中提取可能需要一段时间。

FFTW 分配例程版本：

void get_PDF_CDF(float rbin,std::string filtertype,std::vector<float> field,std::vector<float>& PDF_x,std::vector<int>& PDF_counts,std::vector<float>& PDF,std::vector<float>& CDF){

/*===========================================================
Here we go to Fourier space for applying smoothing filters
=============================================================*/
#ifdef smoothing_filters

fftwf_complex * Fourierfield;
float * field_to_fourier;
Fourierfield = (fftwf_complex*) fftwf_malloc(pow(NGRID,2)*(NGRID/2 + 1)*sizeof(fftwf_complex));
field_to_fourier = (float*) fftwf_malloc(pow(NGRID,3)*sizeof(float));

std::cout << "The FFTW multi-threaded calculations are " << (fftw_init_threads() ? "working" : "NOT working")<< std::endl;

fftw_plan_with_nthreads(nthreads);
fftwf_plan plan_forward;
fftwf_plan plan_backward;
std::cout << "Creating plans for forward and backward Fast Fourier transform..." << std::endl;
plan_forward=fftwf_plan_dft_r2c_3d(NGRID,NGRID,field_to_fourier,Fourierfield,FFTW_ESTIMATE);
plan_backward=fftwf_plan_dft_c2r_3d(NGRID,FFTW_ESTIMATE);

std::copy(field.begin(),field.end(),field_to_fourier);
std::cout << "Executing the forward Fourier transform..." << std::endl;
fftwf_execute(plan_forward);
/*========================
Applying the filter here
========================*/
std::cout << "Applying the " << filtertype << " filter... " << std::endl;
if(filtertype.compare("Gaussian")==0){
        gaussian_f(NGRID,fsize,Fourierfield);
        }    
else if(filtertype.compare("TopHat")==0){
        top_hat_f(NGRID,Fourierfield);
        }

std::cout << "Executing the backward Fourier transform..." << std::endl;
fftwf_execute(plan_backward);
fftw_plan_with_nthreads(1);
fftwf_destroy_plan(plan_forward);
fftwf_destroy_plan(plan_backward);
fftw_cleanup_threads();

std::copy(field_to_fourier,field_to_fourier+FFTnorM,field.begin());

fftwf_free(Fourierfield);
fftwf_free(field_to_fourier);

//RESCALE THE FIELDS (in 3D FFT transform,a factor of NGRID^3 is gained upon returning to the real space)
std::for_each(field.begin(),[&](float& x){ x/=float(FFTnorM);} );
#endif //endif for smoothing filters...
}

std::vector 版本：

void get_PDF_CDF(float rbin,std::vector<float>& CDF){

std::cout << "There are currently " << nthreads << " threads working in parellel" << std::endl;

/*===========================================================
Here we go to Fourier space for applying smoothing filters
=============================================================*/
#ifdef smoothing_filters

std::vector<fftwf_complex> Fourierfield(pow(NGRID,2)*(NGRID/2 +1));

std::cout << "The FFTW multi-threaded calculations are " << (fftw_init_threads() ? "working" : "NOT working")<< std::endl;

fftw_plan_with_nthreads(nthreads);
fftwf_plan plan_forward;
fftwf_plan plan_backward;

std::cout << "Creating plans for forward and backward Fast Fourier transform..." << std::endl;

plan_forward=fftwf_plan_dft_r2c_3d(NGRID,field.data(),Fourierfield.data(),FFTW_ESTIMATE);

std::cout << "Executing the forward Fourier transform..." << std::endl;
fftw_plan_with_nthreads(1);
fftwf_execute(plan_forward);

/*========================
Applying the filter here
========================*/
std::cout << "Applying the " << filtertype << " filter... " << std::endl;

if(filtertype.compare("Gaussian")==0){
        gaussian_f(NGRID,Fourierfield);
        }

else if(filtertype.compare("TopHat")==0){
        top_hat_f(NGRID,Fourierfield);
        }

std::cout << "Executing the backward Fourier transform..." << std::endl;
fftwf_execute(plan_backward);
fftwf_destroy_plan(plan_forward);
fftwf_destroy_plan(plan_backward);
fftw_cleanup_threads();


//RESCALE THE FIELDS (in 3D FFT transform,[&](float& x){ x/=float(FFTnorM);} );

#endif //endif for smoothing filters...
}

过滤器代码如下所示：

void gaussian_f(int NGRID,float fsize,fftwf_complex* Fourierfield){


int i,j,k;
float x_i,y_j,z_k,krnorm,freq,temp =0.0;

freq = fsize*2.0*M_PI/(float)NGRID;

/*=======================================
The input data is real,therefore its
Fourier transform is symmetric
========================================*/

#pragma omp parallel for collapse(3) private(z_k,x_i,temp) num_threads(nthreads)
for( k=0 ; k < int(NGRID/2)+1 ; k++){
        for( j=0 ; j < NGRID ; j++ ){
                for( i=0 ; i < NGRID ; i++ ){

                        z_k = (float)k - 0.5;
                        y_j = (float)j - 0.5;
                        x_i = (float)i - 0.5;
                        krnorm = (x_i*x_i+ y_j*y_j+ z_k*z_k) * freq * freq;
                        temp = exp(-0.5*krnorm);

                        Fourierfield[LOC(i,k)][0] *= temp; Fourierfield[LOC(i,k)][1] *= temp;
       }
     }
   }

}

哪里

#define LOC(ix,iy,iz)  (((iz)*NGRID*NGRID + (iy)*NGRID + (ix))    )

编辑：

我现在绝对确定应用过滤器时会发生错误。将代码的这些部分注释掉允许代码无限期地运行。目前我不知道是什么导致了这个问题。这些函数在程序的顶部声明。现在，我使用的是 gaussian_f，但是 top_hat_f 也会出现同样的问题。此外，现在从 1 个线程执行和删除计划。

void top_hat_f(int NGRID,fftwf_complex* Fourierfield);
void gaussian_f(int NGRID,fftwf_complex* Fourierfield);

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）

allocation c++fftw segmentation-fault