How to create 3D matrix in Eigen

[JamieAl]

I am trying to build a 3D matrix in C++ using Eigen. I read online that instead of 3D matrix, I can create a Vector of 2D matrices. Which will help with performance, especially that I will be using large matrices. So, this is what I attempted, but it's not really working. I think I am misunderstanding something:

static const int nx = 10;
static const int ny = 10; 
static const int nz = 10;

Eigen::MatrixXd W((ny+1),nx); 
W.setZero();
Eigen::MatrixXd M2((ny+1),nx); 
M2.setZero();
vector<Matrix<double, Dynamic, Dynamic>> M3;
M3 W(nx, M2(ny, nz)); // matrix with dimension (k,m,n) or ((nz+1),ny,nx)) 

I mean this doesn't work really, it seems like the wrong syntax, but I would like to know how to correctly use a vector to build 2D matrices and then maybe use a for loop to build a 3d matrix?

Note: I read about this unsupported feature "tensor" in Eigen, but I don't know if using something that's still unsupported is a good idea.

[jonnin]

you want this, maybe:
Eigen::MatrixXd W((ny+1),nx);
vector<Eigen::MatrixXd> ThreeD; //this may need constructor parameters. It may be easier to make it a vector of pointers so you can swap slices faster or split it up for threading etc, depends on what you need here...

for(..)
{
fill in each location of threed with a 2d Eigen::MatrixXd W((ny+1),nx);
}

however I don't think eigen has any functions for 3-d so the created object may be just as well off to be a 3d vector of doubles or possibly better, 1D morphed as 3-d. Are you going to use anything from eigen on 2d slices of it?

[JamieAl]

@jonnin
Thanks! So, I have been reading on the Eigen library tensor and it's sort of nice to use. For example I can build a 3D matrix (tensor of rank 3):

Eigen::Tensor<double, 3> ttt(4,2,2); 
ttt.setZero();
//slicing test: access first matrix in tensor: 
std::array<long,3> offset = {0,0,0};         //Starting point:
std::array<long,3> extent = {1,2,2};  //Finish point:(row,column,matrix)
std::array<long,2> shape2 = {2,(2)};  
std::cout <<  ttt.slice(offset, extent).reshape(shape2) << std::endl;  //Extract slice and reshape it into a 10x10 matrix.
	

This works for me, but I am so new to the syntax that I am struggling to like extract a slice and reshape it to be 4-by-2 matrix (rows=4, cols=2) from the code above. I also, have no real idea of whether using tensor is better than something like this using Eigen:
https://studywolf.wordpress.com/2012/09/16/n-dimensional-matrices-in-c/

[jonnin]

I do not know about tensor and its slice function; you will have to dig on the web for examples.
I will warn you that 'tensor' is a thing in linear algebra and your library may be geared toward this specialty and not general purpose 3d matrix work. You may find it is or is not what you really needed.

if its the same thing, here is their example on slice:

Eigen::Tensor<int, 2> a(4, 3);
a.setValues({{0, 100, 200}, {300, 400, 500},
             {600, 700, 800}, {900, 1000, 1100}});
Eigen::array<Eigen::Index, 2> offsets = {1, 0};
Eigen::array<Eigen::Index, 2> extents = {2, 2};
Eigen::Tensor<int, 2> slice = a.slice(offsets, extents);
cout << "a" << endl << a << endl;
=>
a
   0   100   200
 300   400   500
 600   700   800
 900  1000  1100
cout << "slice" << endl << slice << endl;
=>
slice
 300   400
 600   700

it says it just chops out the region you defined, so 3d to 2d would be something like x1 to x2, y1 to y2, z=constant
or all the rows, all the columns, and a specific depth. You will have to tinker up the syntax to make it do that. Start by verifying you can slice off what you want, then try to reshape from there, don't do it all at once.

[lastchance]

@JamieAl,
Eigen is primarily designed to do linear algebra well. For that reason its principal objects are 2-d matrices and 1-d vectors (in the mathematical sense).

If you want 3-d objects (say, to simulate xyz space) then it is not a great idea to use vectors of matrices, because one dimension (say, z) is being treated very differently from the other two (say, x and y). In 3-d space that is very unnatural - there is nothing sacrosanct about one direction over the others.

Why exactly do you need 3-d objects? Why do you need to "slice" and "reshape" them?

Your end use should guide the choice of tools to be used. I'm not sure that you have ever explained the purpose of the algebraic tools that you have picked.

[JamieAl]

Thanks everyone!
I am trying to run 3D simulations. I have succeeded to write a code in MATLAB that I am trying to convert to C++. The reason that I am trying to slice the 3D matrix/tensor is I wanted to be in a similar shape in MATLAB for easier and better comparisons. So for example in MATLAB I might have something like:

Nx = 4;  
Ny = 4;  
Nz = 4; 

x   = (0:Nx-1)/Nx*2*pi;    % x coordinate in Fourier, equally spaced grid
y   = (0:Ny-1)/Ny*2*pi;
zgl = -cos(pi*(0:Nz)/Nz)'; %Gauss-Lobatto chebyshev points
 
%make mesh
[X,Z,Y] = meshgrid(x,zgl,y);

Lx = 2*pi; 
Ly = 2*pi; 
dx = Lx/Nx;
dy = Ly/Ny;
A = 2*pi / Lx;
B = 2*pi / Ly;

u = Z.^2 .* sin(A*X) .* sin(B*Y); 

Then the size of u is 5-by-4-by-4 and u can be expressed in slices as val(:,:,1) with size 5-by-4. So, I wanted to acess the first matrix similarly, but I am failing to do this val(:,:,1) in C++ using tensors.

[JamieAl]

@jonnin
Thanks! The example you provided is VERY helpful. Although I am really failing to extract say a single column or a column vector for a test. I tried:

Eigen::Tensor<int, 2> a(4, 3);
a.setValues({{0, 100, 200}, {300, 400, 500},
             {600, 700, 800}, {900, 1000, 1100}});
Eigen::array<Eigen::Index, 2> offsets = {0, 0};
Eigen::array<Eigen::Index, 2> extents = {3, 0};
Eigen::Tensor<int, 2> slice = a.slice(offsets, extents);
cout << "a" << endl << a << endl;

which returns nothing for the slice output. I don't understand maybe I am confused. But I thought {0, 0}; and {3, 0}; means first row, column and 4th row and first column which should return a 1-by-4 column vector??

[lastchance]

One of your extents is 0. No size: no output.

[JamieAl]

@lastchance
oh, clearly I have gone mad! Thanks.

[jonnin]

oh, clearly I have gone mad!

Don't worry, you will get used to it, and are in good company.

[JamieAl]

@jonnin
Thanks again.