Incorrect results: mapping c++ array to Eigen Matrix

July 1, 2022

I have used the map feature before to map existing memory into Eigen matrices, however when trying to map an fftw c++ array, I am getting wrong results, almost if a portion (a slice?) of the array is being mapped into an Eigen matrix and not the entire memory block. This is the code I am using:

static const int nx = 10;
static const int ny = 10; 
static const int nyk = ny/2 + 1;
static const int nxk = nx/2 + 1;
static const int ncomp = 2;

fftw_complex *uhk; // this is pert Te
uhk= (fftw_complex*) fftw_malloc((((nx)*(ny+1))*nyk)* sizeof(fftw_complex)); 
memset(uhk, 42, (((nx))*nyk)* sizeof(fftw_complex));


    for (int i = 0; i < nx; i++){
        for (int j = 0; j < nyk; j++){
            for (int k = 0; k < ncomp; k++){
                uhk[i + nyk*j][k] = //taking fft of some expression
        }

    }
    

Eigen::Map<Eigen::MatrixXcd, Eigen::Unaligned> uhOut(reinterpret_cast<std::complex<double>*>(uhkOut),nyk,nx);

 std::cout << uhOut<< '\n';

The results I am getting are,

 (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)
 (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.58,-0.28) (0.95,-0.46) (0.95,-0.46) (0.58,-0.28) (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.59,-0.09) (0.95,-0.14) (0.95,-0.14) (0.59,-0.09) (0.00,-0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)

But, the correct result should be:

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

Am I misunderstanding the point of map here? Why are my results sliced like that?

>Solution :

You’re computing your offsets wrong. [i+nyk*j] is not correct. The largest value of j is (nyk-1). I think you meant for it to be [i*njk+j].

To diagnose this sort of problem, I suggest doing

uhk[j + nyk*i][k] = (k==0)?i:j; // Note I swapped i/j in the offsets

This will create a matrix where the real portions match the column number, and the imaginary portions are the row number. A simple pattern like this makes it obvious the values aren’t being correctly written the expected memory locations.