1D FFT of 3D Eigen tensor along a specific direction

Question

I have the following code in MATLAB where I am trying to reproduce in C++ using the FFTW library:

MATLAB code:

Nx = 8;  
Ny = 6; 
Nz= 8;
Lx =6;
Ly = 6;
xi_x =  (2*pi)/Lx;
yi_y =  (2*pi)/Ly;
xi  = ((0:Nx-1)/Nx)*(2*pi);
yi  = ((0:Ny-1)/Ny)*(2*pi);
x =  xi/xi_x;
y =  yi/yi_y;

zlow = 0;%a
zupp =6; %b
Lz = (zupp-zlow);
zgl = (1/2)*(((zupp-zlow)*zgl) + (zupp+zlow));
[X,Z,Y] = meshgrid(x,zgl,y);
%ICs
A = 2*pi / Lx;
B = 2*pi / Ly;
u = (Z-zlow) .* (Z-zupp) .* sin(A*X).* sin(B*Y);
%Take FFT along only X: 
uh =fft(u,[],2);
%Take FFT along only Y direction:
uh =fft(u,[],3); %since meshgrid(x,z,y) then #3 is Y direction

I am able to get the correct results in C++ for the X direction only so:

uh =fft(u,[],2);

which is:

fftw_plan const the_plan = fftw_plan_many_dft_r2c(
1,  //rank-one transforms
&nx, 
(nz+1), //one transform per row
Tensor2D.data(), //extract 2D matrix out of 3D tensor
nullptr,
1,
reinterpret_cast<fftw_complex*>(Tensor2DK.data()), //output
nullptr,
(nz+1),
1,
FFTW_MEASURE);

Where I define a 3D Tensor and extract a 2D matrix as the following:

static const int nx = 8; 
static const int ny = 6;  
static const int nz = 8; 

Eigen::Tensor<double, 3> uFun((nz+1),nx,ny);   //tensor(row,col,matrix)
uFun.setZero();

Eigen::Tensor<double, 2> Tensor2D((nz+1),nx,ny); 
Eigen::Tensor<double, 2> Tensor2DK((nz+1),nx,ny); 

Tensor2D = uFun.slice(std::array{0,0,1}, std::array{nz+1,nx,1}.reshape(std::array{nz+1,nx});

//Take FFT
fftw_excute(the_plan);

Above code returns correct results of uh =fft(u,[],2); However, I am really unable to take the 1D FFT along only Y similarly as uh =fft(u,[],3); I keep getting results that don't match MATLAB. How can I fix this?

EDIT:

I have simplified my problem more by testing a 2D matrix first since my main issues is I am unable to take the 1D FFT of a 2D matrix along each row using FFTW library. I have a simple MATLAB code that I am using to benchmark my results and it looks like:

Nx=8;
Ny=8;
u = eye(Ny+1,Nx); %very simple choice for now, returns ny+1 rows by nx columns
%take FFT along columns
urow = fft(u,[],1);
ucol = fft(u,[],2); 

Now, in C++ I have tried the following:

1. Take 1D FFT of column-major 2D matrix along each column (this works)

#include "fftw3.h"
#include<unsupported/Eigen/CXX11/Tensor>

static const int nx=8;
static const int ny=8;

using namespace std;
using namespace Eigen;
int main{

Eigen::MatrixXd u(ny+1,nx);
u.setIdentity();
Eigen::MatrixXcd uOut(ny+1,nx);
//1D FFT along cols
fftw_plan const fwd = fftw_plan_many_dft_r2c(
1,
&nx,
u.data(),
nullptr,
(ny+1),
1,
reinterpert_cast<fftw_complex*>(uOut.data()),
nullptr,
(ny+1),
1,
FFTW_MEASURE);
fftw_excute(fwd);
std::cout<<uOut<<endl; //this is correct
}

The above code is correct and returns a result that matches MATLAB: ucol = fft(u,[],2); and output is:

ucol =

   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i
   1.0000 + 0.0000i   0.7071 - 0.7071i   0.0000 - 1.0000i  -0.7071 - 0.7071i  -1.0000 + 0.0000i  -0.7071 + 0.7071i   0.0000 + 1.0000i   0.7071 + 0.7071i
   1.0000 + 0.0000i   0.0000 - 1.0000i  -1.0000 + 0.0000i   0.0000 + 1.0000i   1.0000 + 0.0000i   0.0000 - 1.0000i  -1.0000 + 0.0000i   0.0000 + 1.0000i
   1.0000 + 0.0000i  -0.7071 - 0.7071i   0.0000 + 1.0000i   0.7071 - 0.7071i  -1.0000 + 0.0000i   0.7071 + 0.7071i   0.0000 - 1.0000i  -0.7071 + 0.7071i
   1.0000 + 0.0000i  -1.0000 + 0.0000i   1.0000 + 0.0000i  -1.0000 + 0.0000i   1.0000 + 0.0000i  -1.0000 + 0.0000i   1.0000 + 0.0000i  -1.0000 + 0.0000i
   1.0000 + 0.0000i  -0.7071 + 0.7071i   0.0000 - 1.0000i   0.7071 + 0.7071i  -1.0000 + 0.0000i   0.7071 - 0.7071i   0.0000 + 1.0000i  -0.7071 - 0.7071i
   1.0000 + 0.0000i   0.0000 + 1.0000i  -1.0000 + 0.0000i   0.0000 - 1.0000i   1.0000 + 0.0000i   0.0000 + 1.0000i  -1.0000 + 0.0000i   0.0000 - 1.0000i
   1.0000 + 0.0000i   0.7071 + 0.7071i   0.0000 + 1.0000i  -0.7071 + 0.7071i  -1.0000 + 0.0000i  -0.7071 - 0.7071i   0.0000 - 1.0000i   0.7071 - 0.7071i
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i

2. Take 1D FFT of column-major 2D matrix along each row (INCORRECT) following with this post: How do I use fftw_plan_many_dft on a transposed array of data? where they said: // ... if data is column-major, set istride=howmany, idist=1 // if data is row-major, set istride=1, idist=N

Eigen::MatrixXd u(ny+1,nx);
u.setIdentity();
Eigen::MatrixXcd uOut(ny+1,nx);
//1D FFT along rows by:
// if data is row-major, set istride=1, idist=N
fftw_plan const fwd = fftw_plan_many_dft_r2c(
1,
&nx,
(ny+1),
u.data(),
nullptr,
1, //istride=1
nx, //idist=N
reinterpert_cast<fftw_complex*>(uOut.data()),
nullptr,
(ny+1),
1,
FFTW_MEASURE);
fftw_excute(fwd);
std::cout<<uOut<<endl; //this is incorrect
}

Above code does NOT work and return only zeros.

3. Convert original column-major 2D matrix to a row-major matrix instead and repeat (1) (also does NOT work)

Eigen::MatrixXd u(ny+1,nx);
u.setIdentity();
urow = u.transpose(); //row-major 2D matrix 
Eigen::MatrixXcd uOutrow(ny+1,nx);
fftw_plan const fwd = fftw_plan_many_dft_r2c(
1,
&nx,
urow.data(),
nullptr,
(ny+1),
1,
reinterpert_cast<fftw_complex*>(uOutrow.data()),
nullptr,
(ny+1),
1,
FFTW_MEASURE);
fftw_excute(fwd);
std::cout<<uOut<<endl; //this is incorrect
}

Above code does NOT work either. I am SO confused with what seems to be a simple task. The FFTW library documentation is so difficult to understand with no examples and so I am left even more confused. I think I would like to focus on my edit and get help with the 2D matrix only since it'd be straightforward to implement with tensors once I get it working.

What am I doing wrong? How can I take a simple 1D FFT along each row such as urow=fft(u,[],1); and the correct output should look like:

urow =

   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i
   1.0000 + 0.0000i   0.7660 - 0.6428i   0.1736 - 0.9848i  -0.5000 - 0.8660i  -0.9397 - 0.3420i  -0.9397 + 0.3420i  -0.5000 + 0.8660i   0.1736 + 0.9848i
   1.0000 + 0.0000i   0.1736 - 0.9848i  -0.9397 - 0.3420i  -0.5000 + 0.8660i   0.7660 + 0.6428i   0.7660 - 0.6428i  -0.5000 - 0.8660i  -0.9397 + 0.3420i
   1.0000 + 0.0000i  -0.5000 - 0.8660i  -0.5000 + 0.8660i   1.0000 + 0.0000i  -0.5000 - 0.8660i  -0.5000 + 0.8660i   1.0000 + 0.0000i  -0.5000 - 0.8660i
   1.0000 + 0.0000i  -0.9397 - 0.3420i   0.7660 + 0.6428i  -0.5000 - 0.8660i   0.1736 + 0.9848i   0.1736 - 0.9848i  -0.5000 + 0.8660i   0.7660 - 0.6428i
   1.0000 + 0.0000i  -0.9397 + 0.3420i   0.7660 - 0.6428i  -0.5000 + 0.8660i   0.1736 - 0.9848i   0.1736 + 0.9848i  -0.5000 - 0.8660i   0.7660 + 0.6428i
   1.0000 + 0.0000i  -0.5000 + 0.8660i  -0.5000 - 0.8660i   1.0000 + 0.0000i  -0.5000 + 0.8660i  -0.5000 - 0.8660i   1.0000 + 0.0000i  -0.5000 + 0.8660i
   1.0000 + 0.0000i   0.1736 + 0.9848i  -0.9397 + 0.3420i  -0.5000 - 0.8660i   0.7660 - 0.6428i   0.7660 + 0.6428i  -0.5000 + 0.8660i  -0.9397 - 0.3420i
   1.0000 + 0.0000i   0.7660 + 0.6428i   0.1736 + 0.9848i  -0.5000 + 0.8660i  -0.9397 + 0.3420i  -0.9397 - 0.3420i  -0.5000 - 0.8660i   0.1736 - 0.9848i