## N-dimensional Matrices in C++

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];
for (int d1 = 0; d1 < sizes; d1++) {
matrix[d1] = new MatrixXf* [sizes];
for (int d2 = 0; d2 < sizes; d2++)
matrix[d1][d2] = new (MatrixXf)(MatrixXf::Ones(sizes,sizes).array() * value);
}
return matrix;
}
```

and call it with:

```int sizes = {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,element);
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; d1++) {
for (int d2 = 0; d2 < sizes; 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.

Tagged , ,

## 6 thoughts on “N-dimensional Matrices in C++”

1. Aaron Schurger says:

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

• travisdewolf says:

Hi Aaron,

thanks for the kind words! Really glad that you found this useful! 😀

• Aaron Schurger says:

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 );

}

• travisdewolf says:

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!!

2. Aaron Schurger says:

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!

3. Aaron Schurger says:

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!