C++

Forum

 
Error: Using EigenFFT in for loop

JamieAl (236)
I am testing EigenFFT with some Eigen tensors and matrices. I have been so far successful with implementing 1D FFTs along different directions and comparing my results to MATLAB. So, for example, my MATLAB code looks like:
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
 
clearvars; clc; close all;
%3D FFT test
Nx = 8;  
Ny = 4;%8; 
Nz= 6;%8;
Lx =16;%2*pi;%
Ly = 6;% 2*pi;%
dx = Lx/Nx; % Need to calculate dx
dy = Ly/Ny;

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 = cos ( pi * ( 0 : Nz ) / Nz )';
zgl = (1/2)*(((zupp-zlow)*zgl) + (zupp+zlow));
[X,Z,Y] = meshgrid(x,zgl,y); %this gives 3d grid with z-by-x-by-y size
%ICs
A = 2*pi / Lx;
B = 2*pi / Ly;
uFun = (Z-zlow) .* (Z-zupp) .* sin(A*X).* sin(B*Y);
%Take 1D FFTs along each direction
uh1 =(fft(u,[],1));%each columns
numRows=Nz+1;numCols=Nx;
   for j= 1:numCols
      uk1(:,j,:) = fft(u(:,j,:));
   end
uh2 =(fft(u,[],2));%each row
   for j= 1:numRows
      uk2(j,:,:) = fft(u(j,:,:));
   end
uh3 =(fft(u,[],3));%over 3rd dimension (i.e. Y)
for i= 1:numRows
   for j= 1:numCols
      uk(i,j,:) = fft(u(i,j,:));
   end
end


Now, in C/C++ I am able to reproduce both uh1 =(fft(u,[],1)); and uh2 =(fft(u,[],2)); really the third dimension that I am struggling with the syntax of Eigen and getting errors.

My C/C++ code looks like:
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
 

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

int main(){
Eigen::Tensor<double, 3> uFun((nz+1),nx,ny);//tensor(row,col,matrix)
uFun.setZero();
for(int x = 0; x< nx; x++){
		    for(int y = 0; y< ny; y++){
                for(int z = 0; z< nz+1; z++){
                uFun(z,x,y) = //define function u same as MATLAB
            }
        }
    }
//define matrix and copy the tensor into it
//map tensor to matrix

Eigen::MatrixXd dummy((nz+1),nx); //to match tensor slices
dummy.setZero();
Eigen::Tensor<double, 2> Tensor2D((nz+1),nx);   //tensor(row,col,matrix)
Tensor2D.setZero();

Tensor2D = uFun.slice(std::array{0, 0, 3}, std::array{nz + 1, nx, 1}).reshape(std::array{nz + 1, nx}) ;
dummy = static_cast<MatrixType<double>>(Eigen::Map<const MatrixType<double>>(Tensor2D.data(), (nz+1), nx));
Eigen::MatrixXcd tempOut1((nz+1),nx);
//uh1= fft(u,[],1);
	for (int k=0; k< dummy.cols(); k++){
		tempOut1.col(k) = fft.fwd(dummy.col(k));		
	}
//uh2= fft(u,[],2);
Eigen::MatrixXcd tempOut((nz+1),nx);
Eigen::VectorXcd row(nx);
	for (int k=0; k< (dummy.rows()); k++){
		fft.fwd(row, dummy.row(k));
		tempOut.row(k) = row;
	}

}

The above code works, but now when trying to reproduce this part from MATLAB code:
1
2
3
4
5
6
 
for i= 1:numRows
   for j= 1:numCols
      uk(i,j,:) = fft(u(i,j,:));
   end
end

I get the error:
1
2
3
 
36)  (14.790,-7.123)  (10.458,-5.036)
test.out: /mnt/c/Users/Documents/eigen-3.4.0/eigen-3.4.0/Eigen/src/Core/Block.h:120: Eigen::Block<XprType, BlockRows, BlockCols, InnerPanel>::Block(XprType&, Eigen::Index) [with XprType = Eigen::Matrix<double, -1, -1>; int BlockRows = 1; int BlockCols = -1; bool InnerPanel = false; Eigen::Index = long int]: Assertion `(i>=0) && ( ((BlockRows==1) && (BlockCols==XprType::ColsAtCompileTime) && i<xpr.rows()) ||((BlockRows==XprType::RowsAtCompileTime) && (BlockCols==1) && i<xpr.cols()))' failed.
Aborted (core dumped)


My attempt in C/C++ looks like:
1
2
3
4
5
6
7
8
9
10
11
 
for (int i = 0; i < dummy.rows(); ++i) {
        for (int j = 0; j < dummy.cols(); ++j) {
            int index =  i * (dummy.cols()) + j; // Flattened index
	//int index =  i * (dummy.rows()) + j; // Flattened index
            VectorXcd slice = dummy.row(index);
            tempOut.row(index) = fft.fwd(slice);

        }
    }

std::cout<<tempOut(index)<<endl;

How can I loop over # of rows and columns for each matrix similarly to the MATLAB code?
Last edited on
Topic archived. No new replies allowed.