C++ array of Eigen dynamically sized matrix

Question

In my application, I have a one-dimensional grid and for each grid point there is a matrix (equally sized and quadratic). For each matrix, a certain update procedure must be performed. At the moment, I define a type

typedef Eigen::Matrix<double, N, N> my_matrix_t;

and allocate the matrices for all grid points using

my_matrix_t *matrices = new my_matrix_t[num_gridpoints];

Now I would like to address matrices whose sizes are only known at run time (but still quadratic), i.e.,

typedef Eigen::Matrix<double, Dynamic, Dynamic> my_matrix_t;

The allocation procedure remains the same and the code seems to work. However, I assume that the array "matrices" contains only the pointers to the each individual matrix storage and the overall performance will degrade as the memory has to be collected from random places before the operation on each matrix can be carried out.

Q0: Contiguous Memory?

Is the assumption correct that

new[] will only store the pointers and the matrix data is stored anywhere on the head?
it is beneficial to have a contiguous memory region for such problems?

Q1: new[] or std::vector?

Using a std::vector was suggested in the comments. Does this make any difference? Advantages/drawbacks of both solutions?

Q2: Overloading new[]?

I think by overloading the operator new[] in the Eigen::Matrix class (or one of its bases) such an allocation could be achieved. Is this a good idea?

Q3: Alternative ways?

As an alternative, I could think of using a large Eigen::Matrix. Can anyone share their experience here? Do you have other suggestions for me?

If you want an array of something with a run time size, use a std::vector. — NathanOliver
The problem is not the run time size of the array, but the run time size of the MatrixXd. I need to provide new[] or std::vector somehow with the size of the matrix elements, which are not known at compile time but known at the time when I create the array. — carlosvalderrama
You can do that with a vector. When you create the vector it has no elements in it. Then you can add elements at run time and specify what size they should be when doing so. I'm not familiar with eigen so I don't know if you can overwrite the default construct MatrixXd in an array with one of a different size. If you can, then an array would work as well. — NathanOliver
@carlosvarerrama: you are confusing the memory allocation for the array of matrices (which you can do in your code, preferably with std::vector) and the dynamic memory allocation done internally by each matrix object in the array. You cannot control how each matrix is going to allocate its internal storage, and you can't preallocate this storage yourself. — PlinyTheElder
@carlosvalderrama: sorry, I don't understand what you mean. To change the way the memory is allocated inside the MatrixXd class, you would need to edit the MatrixXd source code (I really don't recommend that you do that). You can allocate an array of matrices, but when you do, each of the matrices in the array will perform its own heap allocation which you can't control if you are just a user of the MatrixXd class. — PlinyTheElder

carlosvalderrama carlosvalderrama · Accepted Answer · 2019-07-26T09:58:21

Let us sum up what we have so far based on the comments to the question and the mailing list post here. I would like to encourage everyone to edit and add things.

Q0: Contiguous memory region.

Yes, only the pointers are stored (independent of using new[] or std::vector).
Generally, in HPC applications, continuous memory accesses are beneficial.

Q1: The basic mechanisms are the same.

However, std::vector offers more comfort and takes work off the developer. The latter also reduces mistakes and memory leaks.

Q2: Use std::vector.

Overloading new[] is not recommended as it is difficult to get it right. For example, alignment issues could lead to errors on different machines. In order to guarantee correct behavior on all machines, use

std::vector<my_matrix_t, Eigen::aligned_allocator<my_matrix_t>> storage;

as explained here.

Q3: Use a large Eigen Matrix for the complete grid.

Alternatively, let the Eigen library do the complete allocation directly by using on of its data structures. This guarantees that issues such as alignment and a continuous memory region are addressed properly. The matrix

Eigen::Matrix<double, Dynamic, Dynamic> storage(N, num_grid_points * N);

contains all matrices for the complete grid and can be addressed using

/* i + 1 'th matrix for i in [0, num_gridpoints - 1] */
auto matrix = storage.block(0, i * N, N, N);