Most efficient option for build 3D structures using Eigen matrices

Question

I need a 3D matrix/array structure on my code, and right now I'm relying on Eigen for both my matrices and vectors.

Right now I am creating a 3D structure using new:

MatrixXd* cube= new MatrixXd[60];
for (int i; i<60; i++) cube[i]=MatrixXd(60,60);

and for acessing the values:

double val;
MatrixXd pos;
for (int i; i<60; i++){
    pos=cube[i];
    for (int j; j<60; j++){
        for (int k; k<60; k++){
            val=pos(j,k);
            //...
        }
    }
}

However, right now it is very slow in this part of the code, which makes me beleive that this might not be the most efficient way. Are there any alternatives?

Daile

If you really need a 3D matrix then there isn't a way around using that nested loop. With a 60x60x60 matrix you've got 216,000 values that you're keeping, so it's always going to be relatively slow. What problem are you trying to solve? Maybe there is a way to solve it without a 3D matrix. – Allegedly

I'm calculating a NxN matrix for a certain number of time steps, and I need to store the data for each instant so that I can perform a weighted average at the end. My question was more about the use of new - more precisely, if there aren't any other types of array (maybe in boost libraries) that can be used to stack several MatrixXd. – Daile

Answer 1

While it was not available, when the question was asked, Eigen has been providing a Tensor module for a while now. It is still in an "unsupported" stage (meaning the API may change), but basic functionality should be mostly stable. The documentation is scattered here and here.

Bivens

Answer 2

An alternative is to create a very large chunk of memory ones, and maps Eigen matrices from it:

double* data = new double(60*60 * 60*60*60);

Map<MatrixXd> Mijk(data+60*(60*(60*k)+j)+i), 60, 60);

At this stage you can use Mijk like a MatrixXd object. However, since this not a MatrixXd type, if you want to pass it to a function, your function must either:

• be of the form foo(Map<MatrixXd> mat)
• be a template function: template<typename Der> void foo(const MatrixBase<Der>& mat)
• take a Ref<MatrixXd> object which can handle both Map<> and Matrix<> objects without being a template function and without copies. (doc)

Gossoon

Answer 3

A solution I used is to form a fat matrix containing all the matrices you need stacked.

MatrixXd A(60*60,60);

and then access them with block operations

A0 = A.block<60,60>(0*60,0);
...
A5 = A.block<60,60>(5*60,0);

Thrashing

Why in these examples do ones use examples with equally spaced dimensions?! Why not 30x40x50?! – Synonym