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

1. new[] will only store the pointers and the matrix data is stored anywhere on the head?
2. 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?

Answer 1

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.

1. Yes, only the pointers are stored (independent of using new[] or std::vector).
2. 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);