SoA support in Eigen

Submitted by Christoph Keller

Assigned to Christoph Keller

Link to original bugzilla bug (#98)
Version: 3.0

Description

In order to fully utilize the VectorEngine using Vectors the data needs to be arranged as a Structure of Arrays.

For example the syntax could be like this.

SoAVector3f<LEN> a(...), b(...), c;
Array<> ar = a.dot(b), ar2 = a.norm()
c = a+b

One could also add some helperfunctions to MatrixBase, but using an explicit SoA class is better in the longrun.

Here is what Gael said on the Mailinglist:

first things first, supporting a new vector engine in eigen is
relatively easy. However, you won't be able to right code
"templatized" on the vector engine, and even with a custom library I
doubt that's doable. The choice of the vector engine has to be
compiler options. So, the way to support multiple vector engines at
runtime (I guess this is what you want?), is to have a trivial
function foo(...) selecting the right implementation function e.g.,
foo_sse(...), foo_sse4(...), etc. The foo_* functions are implemented
only once, by putting it in its own file which will be compiled
multiple times with different compiler options (e.g., -msse, -msse4,
etc.). The actual of the function (foo_sse) will be built using the
preprocessor, e.g.:

EIGEN_CAT(foo_,VECTOR_ENGINE_SUFFIX)(...) { .... }

Now regarding the special interleaved packing of the data, there is
currently no such thing in eigen, however, I think you can easily add
that on top of, e.g., a
Array<float,Dynamic,ei_packet_traits<float>::size,RowMajor>. It will
be initialized with dimension *
((size/ei_packet_traits<float>::size)+((size%ei_packet_traits<float>::size)==0
? 0 : 1) rows where dimension is three in your example, and size is
the number of element. You can easily get the i-th element as follow:

underlying_array.block<dimension,1>((i/ei_packet_traits<float>::size)*dimension,
i%ei_packet_traits<float>::size) = Vector3f(x,y,z);

The idea would be to add a class warping this underlying_array to make
it convenient to use. The main questions are what is the set of
features which have to be supported? How do we want to use it? through
a manual loop over the (e.g., 3x4) blocks? through functors? through
high level expression template code? etc.