Hello David,
Thanks for the pointer, I will have a look at it!
Michael
On 07/29/2019 03:49 PM, David Tellenbach wrote:
Hi Michael,
The approach to use Eigen data structures for the complete grid
sounds very promising -- I will try out this one.
If you decide to got his road you could also use the unsupported
Tensor
module(https://eigen.tuxfamily.org/dox-devel/unsupported/eigen_tensors.html).
If you use Eigen::Tensor<Scalar, 3> each grid matrix could be stored
as a tensor slice along the 3. dimension. You can then access grid
matrix i by using tensor.chip(i, 2) which can be used as an lvalue.
This could be easier in terms of abstraction.
However, some caveats: Even though chip removes a dimension its result
is not a matrix type and a lot of algorithms Eigen provides for dense
matrices won't be available. You could probably fix that by using
Eigen::Map constructed from tensor.data(). If you only perform basic
operations on the grid matrices such as multiplication or addition the
tensor module itself would be sufficient.
Finally to address a question from your very first email:
Is having a contiguous memory region actually beneficial or does it
not make any difference?
A cache friendly memory layout is absolutely crucial for good
performance and memory accesses are often a performance bottleneck in
modern applications.
Cheers,
David
On 29. Jul 2019, at 14:51, Michael Riesch <[email protected]
<mailto:[email protected]>> wrote:
Dear Janos,
Thank you for your answer. The approach to use Eigen data structures
for the complete grid sounds very promising -- I will try out this one.
Best regards,
Michael
On 07/24/2019 05:11 PM, Janos Meny wrote:
My idea for second point is wrong, sorry for that. But I think you
could do the following: You could construct a matrix of type
MatrixXd storage(N, k*N);
and then map the storage matrix like so
auto matrix = m.block(0, i*N, N, N); // i + 1 'th matrix for i
smaller than k
Best regards
Janos
Am Mi., 24. Juli 2019 um 16:59 Uhr schrieb Janos Meny
<[email protected] <mailto:[email protected]>>:
Hey Michael,
I think you should be careful writing something like :
typedef Eigen::Matrix<double, N, N> my_matrix_t;
my_matrix_t *matrices = new my_matrix_t[num_gridpoints];
If you are using c++x for some x smaller than 17, the alignment
of the pointer matrices will, I think, not necessarily be
correct. One option is
to use std::vector<my_matrix_t,
Eigen::aligned_allocator<my_matrix_t>>.
Concerning the second point, one idea could be to write your own
stack allocator (see for example this
<https://howardhinnant.github.io/stack_alloc.html>) and then
pass it as a second argument to std::vector.
But I don't think the performance benefits will be worth the hustle.
Best regards
Janos
Am Mi., 24. Juli 2019 um 16:34 Uhr schrieb Michael Riesch
<[email protected] <mailto:[email protected]>>:
Hello all,
In my simulation tool I have a 1D grid and for each grid
point there is
an equally sized and quadratic matrix. For each matrix a
certain update
operation is subsequently carried out. Currently, the matrix
size is
known at compile time, i.e., something like
typedef Eigen::Matrix<double, N, N> my_matrix_t;
Using this type, I allocate the memory 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, i.e.,
typedef Eigen::Matrix<double, Dynamic, Dynamic> my_matrix_t;
(but still quadratic). 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.
Does anyone have experience with such a use case? Is having
a contiguous
memory region actually beneficial or does it not make any
difference?
Can I allocate memory for the complete array of matrices so
that I have
a contiguous memory region and, if yes, how?
Thanks in advance for your comments!
Best regards,
Michael
