I’m working on coding up an algorithm in C++, based on a Matlab implementation I wrote, and I’ve been using the Eigen library (http://eigen.tuxfamily.org/) for my linear algebra. It’s been relatively straightforward to use so far, but unfortunately it doesn’t (as far as I can tell) have any built in support for N-dimensional matrices. In fact, it appears that N-dimensional matrices aren’t really addressed by any of the standard C++ linear algebra libraries. So, to get my 4-dimensional matrix, I have to make an array of pointers pointing to an array of pointers pointing to my Eigen MatrixXf objects. I initialize my 4-D matrices using this function:

MatrixXf*** initialize_4dMatrix(MatrixXf*** matrix, int sizes[], float value=0.0) { matrix = new MatrixXf** [sizes[0]]; for (int d1 = 0; d1 < sizes[0]; d1++) { matrix[d1] = new MatrixXf* [sizes[1]]; for (int d2 = 0; d2 < sizes[1]; d2++) matrix[d1][d2] = new (MatrixXf)(MatrixXf::Ones(sizes[2],sizes[3]).array() * value); } return matrix; }

and call it with:

int sizes[4] = {1,2,3,4}; MatrixXf ***m = initialize_4dMatrix(m, {1,2,3,4});

To then access the matrices, it’s `(*m[d1][d2])(d3,d4);`

which isn’t bad for accessing one matrix at a time, but what if I want the jth element of the ith row of every matrix (i.e. element (i,j))? The Python access method `m[:,:,i,:]`

(or just `m[:,:,i]`

) would be ideal. Unfortunately, I don’t think there’s any way around writing another function to do it:

static VectorXf* get_elements_across_matrix_array(MatrixXf** matrix, int size, int element[]) { VectorXf* vector = new VectorXf(size); for (int i = 0; i < size; i++) (*vector)[i] = (*matrix[i])(element[0],element[1]); return vector; }

And I would need to write another function that returns a MatrixXf pointer if I wanted to access the ith rows of each matrix.

To delete these matrices I also need to pass into my delete function the sizes of the array. It’s a pain that it can’t dynamically discover the array size, but I don’t think there’s any way around it in C++.

static void delete_4dMatrix(MatrixXf*** matrix, int sizes[]) { for (int d1 = 0; d1 < sizes[0]; d1++) { for (int d2 = 0; d2 < sizes[1]; d2++) delete matrix[d1][d2]; } }

Additionally, I need to have a separate function for initializing my 3D matrices, and any other dimension that I would need the way I have it set up now. Less than ideal. My guess is that there’s a elegant way to handle this with void pointers, but at the moment I only need 3D and 4D matrices so it’s not worth the effort for me to set it up.

So there are a few obstacles to using N-dimensional matrices with Eigen, but it’s really not a problem once you get a handle on it and how to properly reference different elements and matrices through your (N-2)-dimensional array of pointers.

Hi Travis,

Just wanted to thank you for this posting, which has been enormously helpful to me as I am just getting to know Eigen and odeint, which I am using to set up and solve ODEs.

Many thanks!

Aaron

Hi Aaron,

thanks for the kind words! Really glad that you found this useful!đŸ˜€

You’re welcome! I’ve made some progress thanks to your posting, now I am just a bit stuck trying to get Eigen and odeint to work together. I am assuming that you can use an Eigen MatrixXd as the “state_type” for odeint (as below). But that doesn’t seem to work. Have you used Eigen and ideint together? Does this look like the proper way to do it?

typedef Eigen::Matrix state_type;

void odefun( const state_type &u , state_type &dudt , double t )

{

…

}

void write_odefun( const state_type &x , const double t )

{

int i;

cout << t << '\n';

for (i=0; i<=9; i++) {

cout << x(i) << '\n';

}

}

int main(int argc, char **argv)

{

…

state_type u(3,3); // initial conditions

srand((unsigned int) time(0));

u.setRandom(); // picks random numbers from -1 to 1

u*=0.1; // make the random numbers smallish

…

integrate( odefun , u , 0.0 , 100.0 , 0.1 , write_odefun );

…

}

Hmm. I’m unfortunately a bit swamped at the moment, and don’t have time to really dive in, but it seems like there may have been similar problems from others: http://stackoverflow.com/questions/22782224/using-boostodeint-with-eigenmatrix-as-state-vector

and

http://stackoverflow.com/questions/28965694/using-a-stdvectoreigenvector3d-as-state-type-in-odeint

if those don’t do it for you I can try to take a look asap but it might be a bit! Please let me know if you get it! Good luck!!

Thanks for the links! The first one seems highly relevant. I’ll mess around with it a bit and see if I can get it working. If I get really desperate, then I may take you up on your offer! I’m not in a huge hurry. Thanks again!

Everything compiles and runs OK, as long as I comment out the critical line…

integrate( dudt_fun , u , 0.0 , 100.0 , 0.1 , write_dudt );

So essentially once I try to use the ode solving capabilities of odeint, I get a compile error:

error: no matching function for call to ‘integrate(void (&)(state_type&, state_type&, double), state_type&, double, double, double, void (&)(const state_type&, double))’

If you have any ideas, I’d be very grateful!