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