FFTW3 segmentation errors while using std::vector container for FFT

  allocation, c++, fftw, segmentation-fault

I would like to ask for your advice regarding the use of std::vector container in performing FFT.

In my program, I use a certain function (get_PDF_CDF) in a loop. When I perform the FFTs using std::vector containers, during a certain iteration (usually the second one), I am getting a segmentation error / memory corruption errors.

I know that FFTW best works with alignment provided by fftw_malloc and fftw_free, but I wanted to spice things up and use std::vector consistently throughout the whole program.

I have rewritten the function to use the recommended allocation routines – it essentially copies the contents of the vector to an array of floats allocated and deallocated using fftw routines – and it works perfectly. It also involves simple rewriting of gaussian_f and top_hat_f filters. However, I would like to learn if I can in any way modify the version with std::vectors so that it does not produce segmentation faults.

Below are the two versions, truncated to their essential content. My apologies if this is not sufficient for error-tracing – in that case I will provide a working example, although it might take a while to extract from the code.

FFTW allocation routines version:

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*) fftw_malloc(pow(NGRID,2)*(NGRID/2 + 1)*sizeof(fftwf_complex));
field_to_fourier = (float*) fftw_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, NGRID, field_to_fourier, Fourierfield, FFTW_ESTIMATE);
plan_backward=fftwf_plan_dft_c2r_3d(NGRID, NGRID, NGRID, Fourierfield, field_to_fourier, 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, fsize, 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();

field.insert(field.end(), &field_to_fourier[0], &field_to_fourier[FFTNORM]);
fftw_free(Fourierfield);
fftw_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(), field.end(), [&](float& x){ x/=float(FFTNORM);} );
#endif //endif for smoothing filters...
}

std::vector version:

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){

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, NGRID, NGRID, field.data(), Fourierfield.data(), FFTW_ESTIMATE);
plan_backward=fftwf_plan_dft_c2r_3d(NGRID, NGRID, NGRID, Fourierfield.data(), field.data(), FFTW_ESTIMATE);

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, fsize, 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();

//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(), field.end(), [&](float& x){ x/=float(FFTNORM);} );

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

Source: Windows Questions C++

LEAVE A COMMENT