How to map/copy complex array of pointers to Eigen matrix?

[JamieAl]

I am able to both copy and map array of pointers to an Eigen vector/matrix and vice versa as the following:

//Example 1:
    std::vector<double> xl(nx+1); 
    for (int i = 0; i< nx; i++){
		xl[i] = (i)*2*EIGEN_PI/nx; 
	}

	Eigen::VectorXd exl = Eigen::Map<Eigen::VectorXd, Eigen::Unaligned>(xl.data(), xl.size());
//Example 2:
Eigen::Matrix< double, (ny+1), (nx)> u; 
	u.setZero(); 
	for(int i = 0; i< nx; i++){
		for(int j = 0; j< ny+1; j++){
			u(j + (ny+1)*i) = pow((eYY(j + (ny+1)*i)-0.5),2) * sin(A * eXX(j + (ny+1)*i)); 
		}
	}
double *uhk;
	uhk = (double*) fftw_malloc((((ny+1)*nx))*sizeof(double)); 
	memset(uhk, 42, (((ny+1)*nx))* sizeof(double));
memcpy(uhk,u.data(),u.size() * sizeof(double));

Both above examples work fine for me, but when I define a complex array of pointers and try to copy it to a complex Eigen matrix it doesn't really work.

fftw_complex *uhkOut; 
	uhkOut = (fftw_complex*) fftw_malloc((((nx)*(ny+1))*nyk)* sizeof(fftw_complex)); 
	memset(uhkOut, 42, (((nx))*nyk)* sizeof(fftw_complex)); 
//initialize uhkout in for loop
//attempts
Eigen::MatrixXcd uh = Eigen::Map<Eigen::MatrixXcd, Eigen::Unaligned>(uhkOut .data(), uhkOut .size()); //ERROR
Eigen::Map<Eigen::MatrixXcd, Eigen::RowMajor> mat(&uhkOut[0], nx, nyk); //ERROR

How can I copy or map this to a complex Eigen matrix correctly?

[mbozzi]

You have to insert a type conversion (a reinterpret_cast conversion). This is okay because of an exception to the rules.

#include <Eigen/Dense>
#include <iostream>
#include <cstdlib>
#include <complex> 

// Just for example's sake:
using my_fftw_complex = double[2];
using my_std_complex = std::complex<double>;

[[nodiscard]] static void* my_fftw_malloc(std::size_t bytes) noexcept
{ if (void* p = std::malloc(bytes)) return p; else std::abort(); } 
static void my_fftw_free(void* p) noexcept { std::free(p); }

[[nodiscard]] static my_fftw_complex* my_fftw_alloc_complex(std::size_t n) noexcept
{ return static_cast<my_fftw_complex*>(my_fftw_malloc(n * sizeof(my_fftw_complex))); }

// See [complex.numbers.general]/4:
[[nodiscard]] static std::complex<double>* to_std_complex(my_fftw_complex* p) noexcept
{ return reinterpret_cast<std::complex<double>*>(p); } 

int constexpr nx = 10;
int constexpr ny = 11;

int main()
{ 
  int const size = nx * ny;

  my_fftw_complex* const buf = my_fftw_alloc_complex(size);

  for (int i = 0; i < size; ++i) 
  { 
    buf[i][0] = i;
    buf[i][1] = 0.0; 
  }

  using stride = Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>; 
  using map = Eigen::Map<Eigen::MatrixXcd, Eigen::Unaligned, stride>;
  map uh(to_std_complex(buf), ny, nx, stride{1, nx});

  uh *= 1.0 / 2.0;

  std::cout << uh << '\n';
}

Live demo:
https://godbolt.org/z/WTzhbqvYj

By the way there is no copying involved: Eigen is just interpreting the memory you've pointed it to as a matrix. This means the memory must not be freed until all matrices that have a map of it are gone.

[JamieAl]

What is stride here mean? also I am not sure what line 40 is doing?
When attempting to follow with this example my complex Eigen matrix doesn't return the same values of uhkOut . Why map is used here differently than my example? Thanks!

[JamieAl]

To shed some light the results I am getting here are:
Output of uhkOut

     (0.00,0.00)     (5.14,0.00)     (8.32,0.00)     (8.32,0.00)     (5.14,0.00)     (0.00,0.00)    (-5.14,0.00)    (-8.32,0.00)    (-8.32,0.00)    (-5.14,0.00)
     (0.00,0.00)    (2.43,-2.35)    (3.93,-3.81)    (3.93,-3.81)    (2.43,-2.35)    (0.00,-0.00)    (-2.43,2.35)    (-3.93,3.81)    (-3.93,3.81)    (-2.43,2.35)
     (0.00,0.00)    (0.52,-1.04)    (0.85,-1.68)    (0.85,-1.68)    (0.52,-1.04)    (0.00,-0.00)    (-0.52,1.04)    (-0.85,1.68)    (-0.85,1.68)    (-0.52,1.04)
     (0.00,0.00)    (0.57,-0.55)    (0.93,-0.89)    (0.93,-0.89)    (0.57,-0.55)    (0.00,-0.00)    (-0.57,0.55)    (-0.93,0.89)    (-0.93,0.89)    (-0.57,0.55)
     (0.00,0.00)    (0.58,-0.28)    (0.95,-0.46)    (0.95,-0.46)    (0.58,-0.28)    (0.00,-0.00)    (-0.58,0.28)    (-0.95,0.46)    (-0.95,0.46)    (-0.58,0.28)
     (0.00,0.00)    (0.59,-0.09)    (0.95,-0.14)    (0.95,-0.14)    (0.59,-0.09)    (0.00,-0.00)    (-0.59,0.09)    (-0.95,0.14)    (-0.95,0.14)    (-0.59,0.09)

Then the Eigen matrix uh

(0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)
(0.58,-0.28) (0.59,-0.09)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00) (0.95,-0.46) (0.95,-0.14)  (0.00,0.00)  (0.00,0.00)
 (0.00,0.00)  (0.00,0.00) (0.95,-0.46) (0.95,-0.14)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00) (0.58,-0.28) (0.59,-0.09)
 (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00) (0.00,-0.00) (0.00,-0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)
 (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)
 (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)

so I see some results are correct but not everything

[JamieAl]

Actually, I see where the problem is, thanks @mbozzi this was significantly helpful.

[mbozzi]

Glad I could help, sorry to reply late.

[JamieAl]

@mbozzi np, I appreciate the help!