1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43
|
void r2cfft3d(Eigen::Tensor<double, 3>& rArr, Eigen::Tensor<std::complex<double>, 3> cArr){
fftw_complex *input_array;
fftw_complex *output_array;
input_array = (fftw_complex*) fftw_malloc(nx*ny*nz * sizeof(fftw_complex));
output_array = (fftw_complex*) fftw_malloc(nx*ny*nz * sizeof(fftw_complex));
for (int i = 0; i < nx; ++i) {
for (int j = 0; j < ny; ++j) {
for (int k = 0; k < nz; ++k) {
{
input_array[k + nz * (j + ny * i)][REAL] = rArr(k,i,j);
input_array[k + nz * (j + ny * i)][IMAG] = 0;
}
}
}
}
fftw_plan forward = fftw_plan_dft_3d(nx, ny, nz, input_array, output_array, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_execute(forward);
fftw_destroy_plan(forward);
fftw_cleanup();
for(int i=0; i < nx; ++i){
for(int j=0; j < ny; ++j) {
for(int k=0; k < nz; ++k) {
cArr(k,i,j).real(output_array[k + nz * (j + ny * i)][REAL]);
cArr(k,i,j).imag(output_array[k + nz * (j + ny * i)][IMAG]);
}
}
}
std::cout<<cArr<<endl; //this is correct
fftw_free(input_array);
fftw_free(output_array);
}
| |