Search code examples
c++tensorfloweigeneigen3

How to copy a c++ array into eigen tensor


I have a code that takes the 3D FFT of an Eigen tensor and returns an Eigen tensor output. I am using a c++ array to pass it to the FFTW plan however I am having an issue getting the correct output (i.e. the size is off). Is there a better way to copy the content of the ''output_array'' into an Eigen tensor using Eigen::map instead of a for loop?

Code example:

Main code

static const int nx = 4;
static const int ny = 4;
static const int nz = 4;
void r2cfft3d(Eigen::Tensor<double, 3>& rArr, Eigen::Tensor<std::complex<double>, 3>& cArr){    
       
        fftw_complex *input_array;
        input_array = (fftw_complex*) fftw_malloc(nx*ny*nz * sizeof(fftw_complex));
        memcpy(input_array, rArr.data(), nx*ny*nz * sizeof(fftw_complex));

        fftw_complex *output_array;
        output_array = (fftw_complex*) fftw_malloc(nx*ny*nz * sizeof(fftw_complex));
        
        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(i,j,k).real(output_array[i + nz * (j + ny * k)][REAL]);
                    cArr(i,j,k).imag(output_array[i + nz * (j + ny * k)][IMAG]);
                
                }

            }
        }
        
        fftw_free(input_array);
        fftw_free(output_array);
    }

This returns the following:

     1e+11        1e+11        1e+11        1e+11        1e+11        1e+11        1e+11        1e+11 -9.53674e-07 -9.53674e-07 -9.53674e-07  2.86102e-06            0            0            0            0
       1e+11        1e+11        1e+11        1e+11        1e+11        5e+11        1e+11        1e+11 -4.76837e-06 -4.76837e-06 -4.76837e-06 -1.62125e-05 -1.52588e-05 -1.52588e-05 -1.52588e-05 -1.52588e-05
       1e+11        1e+11        1e+11        1e+11        1e+11        1e+11        1e+11        1e+11 -9.53674e-07 -9.53674e-07 -9.53674e-07  2.86102e-06            0            0            0            0
       1e+11        1e+11        1e+11        1e+11        1e+11        1e+11        1e+11        1e+11 -8.58307e-06 -8.58307e-06 -8.58307e-06 -4.76837e-06 -7.62939e-06 -7.62939e-06 -7.62939e-06 -7.62939e-06

While the correct answer is

1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11
1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11
1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 5e+11 1e+11 1e+11 1e+11 1e+11 1e+11
1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11 1e+11

It seems like some columns are being filled with zeros for some reason and the rows and columns are off as well (i.e. location of 5e11)

Is there a way to use something like this instead of the for loop here:

TensorMap<Eigen::Tensor<std::complex<double>, 3>> test(output_array, nx,ny,nz);

Solution

  • I was able to fix this by using memcpy actually. So, instead of the for loop I added the line:

    memcpy(cArr.data(),output_array, nx*ny*nz * sizeof(fftw_complex));