Re: [eigen] Initial implementation of tensor support

On Wed, Sep 25, 2013 at 12:32 PM, Christian Seiler <christian@xxxxxxxx> wrote:

Dear all,

I've needed to store objects with more than two indices, i.e. tensors
of arbitrary rank. Unfortunately, Eigen doesn't support those so far.

Since I needed this myself, I decided to write an own tensor class in
C++. But I feel this could also be useful for other people, so I
decided to contribute this to Eigen.

Please note that so far I have only needed storage and no other cool
features Eigen provides for matrices and vectors. For this reason,
there is currently only a single tensor class and no support for
assignments, adding/subtracting tensors, no _expression_ engine etc.
It does require C++11 extensively, so g++ 4.6, clang++ 3.1 or Intel's
compiler 13.1 is needed. (I didn't have a chance to test with MSVC
yet.)

I have adapted my own class a bit to the way Eigen's other code is
structured, but not completely yet. My plan was, since I only have
a limited amount of time, to send my current (working, but incomplete)
version now, and then try to improve upon that later on.

You can find it under
<https://bitbucket.org/chris-se/eigen/src/760259d391f33466e36d19795ee4bc84040dce3f/?at=tensor>
(I'm more used to Git than Mercurial, so please forgive me if I don't
know the proper workflow that you use here.)

Currently, this is added to the unsupported/ section. I've provided
unit tests for nearly everything that's currently implemented. In order
to compile the stuff, you need to pass --std=c++0x for at least g++,
(and probably icpc), so the tests are only activated if you define
cmake -DEIGEN_TEST_TENSOR. They successfully work with g++ 4.6 and 4.7
as well as clang++ 3.2 and icpc 13.1 (i.e. Composer XE 2013.5).

Current features are:

- Create a tensor object with sizes:
Tensor<double, 3> myTensor(42, 23, 10);
Tensor<double, 4, RowMajor> myRowMajorTensor(5, 6, 12, 5);
(This uses Eigen's initialization logic, i.e.
EIGEN_INITIALIZE_COEFFS_IF_THAT_OPTION_IS_ENABLED)

- Copy/Move constructors

- Setting it to zero:
myTensor.setZero();

- data() / size() accessors

- Single-coefficient access
tensor(i) (0 <= i < size())
tensor[i] (only for rank 1 tensors)

- Multi-coefficient access (load & save)
tensor(i,j,k)

- Storage considering symmetries (see below)

Not implemented yet, but I think those things would be really cool and
I will try to implement this one thing at at time at some point in the
future:

- Kind of find a more portable way to make CMake compile the tests
with --std=c++0x for g++... and maybe skip the test if the compiler
has no chance of supporting this
- Packet access
- Assignment of tensors
- Static constructors as for matrices (::Random, ::Zero, ...)
- Assignment operators
- _expression_ engine (A = B + C + D coefficientwise)
- Compile-time size tensors
- Sub-tensors otherTensor = tensor(_1, _2, _3, 42)
(syntax only an example)
- for tensors and/or tensor expressions of rank 1/2 methods to convert
them to matrices (.matrix()) so that they may be used in matrix
expressions (e.g. tensor(_1, _2, 42).matrix().eigenValues())
- Efficient tensor contractions
- possibly more...?

-------------------------------------------------------------------------
Side note: Tensor index symmetries
(aka crazy feature that is already implemented now)

I've already implemented a quite complex feature now (because I needed
it myself): One may specify symmetry relations between tensor indices
and cause the tensor class to set all coefficients according to that
symmetry. For example, you could create the Levi-Civita-Symbol in 3
dimensions just by:

SGroup<3, AntiSymmetry<0,1>,AntiSymmetry<1,2>> s3;
Tensor<int, 3> epsilon(3,3,3);
epsilon.setZero();
epsilon.symCoeff(s3, 0, 1, 2) = 1;

The storage still stores all coefficients, but all the five other
coefficients will now be set to the correct values. Note that this
uses C++ meta templates to actually resolve the symmetry group at
compile time so that the loop to set the other permutations of those
coefficients (including the correct sign) will be unrolled at compile
time, so the above code should really be equivalent to:

Tensor<int, 3> epsilon(3,3,3);
epsilon.setZero();
epsilon.coeff(0, 1, 2) = 1;
epsilon.coeff(1, 2, 0) = 1;
epsilon.coeff(2, 0, 1) = 1;
epsilon.coeff(1, 0, 2) = -1;
epsilon.coeff(0, 2, 1) = -1;
epsilon.coeff(2, 1, 0) = -1;

If the symmetry group is too large (you can cause most compilers to
segfault if you try to use that template implementation with too large
a group) there is some logic that will choose to generate the group at
runtime instead.
-------------------------------------------------------------------------

Anyway: I think a tensor class itself might really be useful for quite
a few people (the symmetry stuff I don't know about, but I decided to
include it anyway) and I do want to continue working on that in the
future (although not in the next month specifically) and wanted to ask
whether you'd consider including this in Eigen.

Christian