I was planning to create a custom Matrix class for my application, but having stumbled upon Eigen and being seriously impressed, I thought to see if there's a way to achieve what I'm trying to do with Eigen.

My application uses some (yet undetermined) type of Genetic Algorithm to select coefficients for many different matrices. My original concept involved the creation of a single, continguous array of (aligned) memory for the parameter vector. Embedded within that large chunk of memory are all of the different matrices used in the computation. I have been intended to utilize Intel MKL and Intel VML for the matrix and array operations, respectively.

A short, simplified illustration (ignoring padding for stride) with small matrices (actual matrices probably between 10k - 10M elements)

Code: Select all

P = [A_11 A_21 A_12 A_22 B_11 B_21 B_31 B_12 B_22 B_32 B_13 B_23 B_33]
A = &P[A_offset]; B = &P[B_offset];


I can then do matrix operations such as

Code: Select all

C = Ax*B ...

I can also perform VML operations on P itself, so that I can have something like P_new = P1 - \alpha * P2.

The above requirement is not overly complicated and a very simple matrix library isn't that hard to implement, but I would really like to see if it's possible to accomplish this with Eigen, to take advantage of the flexibility afforded by the implementation.

I know it is possible to map Eigen to existing arrays; however, these would (in my mind) be allocated manually with Intel MKL custom malloc and I am somewhat concerned above incompatibilities which may come as a result.

Update:
To clarify the issue:
I have 2 different tasks to do with the same data:

Array operations as part of the Genetic Algorithm

Matrix operations which are part of the function(s) being optimized

Of these two, the second should be far more computationally intensive and would benefit most from Eigen (or BLAS). My expectation is that the first task will involve somewhere on the order of 10,000 operations (in the form d = a + b*(c-d) ) on arrays with approximately 1M-2M elements. I know that I could keep a collection of separate matrices as the parameter "vector" and then process the collection in parallel (no race conditions since each matrix is completely independent of the others). However, this does not seem to be graceful, and I don't believe it is the best way to do it. I do believe that something similar to what I described above would be the ideal, but I need help to understand if it can be done with Eigen.


Any insights would be greatly appreciated.
Thanks, Shmuel

What you want is the Map class: http://eigen.tuxfamily.org/dox/group__T ... Class.html, e.g.:

Map<MatrixXd> A(P+A_offset, rows, cols);

Then you can utilize A as an Eigen's matrix (read/write). If you need strides:

Map<MatrixXd, 0, OuterStride<> > A(P+A_offset, rows, cols, OuterStride<>(lda) );

@ggael,
Thank you for your comment. I did stumble on the Map class just prior to my asking the question (or to be a bit more specific, I stumbled on the OuterStride class), which does appear to generally address my question.

Prior to finding Eigen, I had originally intended to write a simple wrapper class to layout the memory as described. I started looking for alternatives after attempting to create a template-based version of my wrapper (for both single/double precision FP), and coming to the realization that my skills with templates in c++ are rather lacking, to say the least...

I believe that the point which leaves me somewhat concerned with the Map class, deals with the memory allocation. Perhaps this is a non-issue, but I do not know enough about Eigen to comment one way or another:

Among Intel's recommendations regarding best practices for performance is the use of their (de-)allocators -- mkl_malloc, mkl_free. I do not know the exact mechanics of these functions and I do not want to have issues due to some incompatibility. I am concerned of a possible situation where Intel's allocator does something unique with the memory layout for the benefit of performance but that is unexpected from the point of view of Eigen.


My inclination is that mkl_malloc simply ensures optimal alignment and that the pointer returned would be no different than from any allocator. If that is the case, then I suspect that Eigen can/will work perfectly. However, I'd really appreciate if you or anyone else has any additional insight into the memory management behaviour of Eigen and possible conflict with that of MKL.

Thanks