C++

Forum

 
How to copy data to an Eigen tensor?

JamieAl (236)
I am trying to copy the output of an array after taking the 3D FFT to a complex Eigen tensor and I am not sure how to do it correctly.

Example:
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
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
 
static const int nx = 4; 
static const int ny = 4;  
static const int nz = 4;

double Lx = 2*EIGEN_PI;
double Ly = 2*EIGEN_PI;
double dx = Lx / nx;
double dy = Ly / ny;

Eigen::Tensor<double, 3> eXX(nx,ny,nz); 
eXX.setZero(); 
Eigen::Tensor<double, 3> eYY(nx,ny,nz); 
eYY.setZero(); 
Eigen::Tensor<double, 3> eZZ(nx,ny,nz); 
eZZ.setZero(); 
for(int i = 0; i< nx; i++){
	for(int j = 0; j< ny; j++){
		for(int k = 0; k< nz; k++){ 
			eXX(k,i,j) = i*dx;
			eYY(j,i,k) = j*dy;
			eZZ(j,i,k) = k*dz; 
    }
  }		
}

double A = (2 * EIGEN_PI)/Lx;
double B= (2 * EIGEN_PI)/ Ly;
Eigen::Tensor<double, 3> uFun(nx,ny,nz); 
uFun.setZero(); 
	
for(int i = 0; i< nx; i++){
   for(int j = 0; j< ny; j++){
      for(int k = 0; k< nz; k++){ 
		uFun(k,i,j) = pow(eZZ(k,i,j),2.0) * sin(A * eXX(k,i,j)) * sin(B * eYY(k,i,j));
			}

		}		
	}

#define IMAG 1
#define REAL 0
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] = uFun(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();

I am trying to copy the output to a tensor with the definition:
1
2
 
Eigen::Tensor<std::complex<double>, 3> Test(nx,ny,nz);  
Test.setZero();
Last edited on
JamieAl (236)
since Test.real() returns the real part of the complex tensor Test. I tried something like this:
1
2
3
4
5
6
7
8
9
10
11
12
 
for(int i=0; i < nx; ++i){

     for(int j=0; j < ny; ++j) {

        for(int k=0; k < nz; ++k) {
            Test(k,i,j).real() = output_array[k + nz * (j + ny * i)][REAL];
            Test(k,i,j).imag() =  output_array[k + nz * (j + ny * i)][IMAG];
               
      }

   }
}

but this give me the error:
expression must be a modifiable lvalue
mbozzi (3932)
Using std::complex::real and imag:
1
2
 
Test(k,i,j).real(output_array[k + nz * (j + ny * i)][REAL]);
Test(k,i,j).imag(output_array[k + nz * (j + ny * i)][IMAG]);
https://en.cppreference.com/w/cpp/numeric/complex/real

Or even 
1
2
3
4
 
[[nodiscard]] constexpr std::complex<double>
to_std_complex(fftw_complex c) { return std::complex<double>{ c[REAL], c[IMAG] }; }
// ...
Test(k, i, j) = to_std_complex(output_array[k + nz * (j + ny * i)]);
Last edited on
JamieAl (236)
@mbozzi
wow! I mean, I almost gave up on this particular issue. Thanks a lot!
Topic archived. No new replies allowed.