Re: [eigen] Multi Dimensional Arrays

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

2011/2/8 Jim Bosch <talljimbo@xxxxxxxxx>:
> On 02/08/2011 05:13 PM, Andy Somogyi wrote:
>>
>> From what I understand, the only way to handle multi dimensional
>> arrays with Eigen is something along the lines of
>>
>> Array<Array<double,Dynamic,Dynamic>  ,Dynamic, Dynamic>  ...
>>
>> This is really far from ideal as it really is an array of arrays
>> (separate memory blocks).  Has there been any talk of merging in some
>> multi dimensional support to Array, say from the now defunct blitz++
>> code, or adding it a new?
>>
>
> That's not entirely true; while Eigen is limited to 2 dimensions,
> Eigen::Matrix is not an array of arrays - they're allocated as a single
> memory block.

But here he's talking about  Array<Array<double,Dynamic,Dynamic>
,Dynamic, Dynamic>  and he's right to say that each
Array<double,Dynamic,Dynamic> would be separately dynamically
allocated.

Using Array as scalar type is only going to be efficient for fixed
sizes (on the inner array type that's used as scalar type). And even
that will only get really efficient in Eigen 3.1, see
http://eigen.tuxfamily.org/bz/show_bug.cgi?id=100.

>
>> I really really really just want to use A SINGLE LIBRARY, and have to
>> rely on Eigen for linear algebra, and boost multi-array or blitz++
>> for multi-dimensional arrays. What would be nice is to have a
>> multi-dim array, be able to grab a dimension or sub-array out of it,
>> and use it directly as an Eigen Array / Matrix.
>>
>
> I feel your pain, and I've been working on a multidimensional array library
> that's designed to let Eigen do all the matrix math:
>
> http://code.google.com/p/ndarray/
>
> While that isn't the complete solution solution you're looking for, at least
> it's a clear division of labor.
>
> Essentially, while you can make an ndarray::Array of any dimension, you can
> view 1- or 2-dimensional subarrays of those as Eigen-conforming objects, or
> construct ndarray views into Eigen::Matrix objects.

This is a good approach, I'd like to encourage you in this direction.
I almost wrote back to the list to suggest that this might be the best
approach to take balancing reasonable performance with reasonable ease
of implementation.

Like with matrices, the hard part is products. The name of the game is
memory locality / cache friendliness. By leveraging Eigen's optimized
matrix product implementation and applying it to 2-dimensional slices
in your N-dim arrays, you should be able to get somewhere. Depending
on the storage order/ cases, you may have to copy slices into
temporaries to have them contiguously stored in memory.

Benoit

>
> That won't get you all the way there if your main application is
> multidimensional tensor contractions, but it does provide a nice interface
> for interfacing with Python/Numpy and moving Eigen-compatible objects around
> with reference counting instead of deep copies.
>
> Good Luck!
>
> Jim Bosch
>
>
>