Boris Jeremic built a neat C++ tensor library called nDarray.
His paper and source code are available online:
http://sokocalo.cee.clarkson.edu/~jeremic/publications/oofepPaper.eps
http://civil.colorado.edu/nDarray/
Here are some initial thoughts about tensors for blitz (and/or Pooma).
1. What is a reasonable upper bound on rank for a tensor?
Blitz++ arrays go up to rank 11; Jim Crotinger pointed out that
a tensor of rank 11 would have at least 2^11=2048 elements, which
makes it unreasonable as a value class. A more reasonable limit
would be rank 4 or 6; 2^6 = 64.
But are there applications in which people would actually
want tensors of rank > 6? Mind you, we're talking *computational*
applications..
2. Many kinds of tensors
There many kinds of tensors; each kind needs different
representation, indexing scheme, and algebraic semantics.
Some examples:
Dense-- your vanilla asymmetric, dense tensor.
Fully symmetric: e.g. A_ijk = A_jki = A_jik. Indexing such
tensors is tricky and involves a lot of combinatorics.
Partly symmetric: the elements are invariant under some permutation
group defined on the indices.
Identity
Kronecker delta
Permutation tensor
Then there are lots of application-specific tensors which have
various oddities (e.g. antisymmetries).
How to handle all these different kinds of tensors? As a single
class with an "engine" template parameter? As a family of classes?
It should be easy for users to create their own tensor
classes, because it's hard to predict all possible varieties
people might need.
3. Should tensor indices be associated with a particular tensor?
In Boris's library, a tensor object remembers its tensor indices.
For example,
tensor I_ijkl = I2("ij")*I2("kl")
The tensor I_ijkl knows that its indices are i,j,k and l. This
is not the case in blitz. Is this important?
4. How to handle contractions
Boris's library does implicit contractions, e.g.
tensor r = A("ik")*B("jk")
results in r_ij = sum over k of A_ik*B_jk
In Blitz contractions are done explicitly, e.g.
Array<float,2> r = sum(A(i,k) * B(j,k), k);
If you're doing lots of tensor notation, Boris's is more compact.
The Blitz style is arguably more flexible, since you can substitute
other operations for "sum".
5. Expression templates?
Steve Karmesin suggested that expression templates aren't necessary,
and that you can just do pairwise evaluation with construct-in-place
via a proxy return type. For simple expressions like
A + B + C
this could be equally or more efficient than expression templates,
particularly on not-so-good compilers. The big problem with
expression templates is that you generate lots of overhead
during parsing; not-so-good compilers won't get rid of this,
so there might be a performance loss if expression templates
were used.
I basically agree with this, but now I'm wondering how to handle
complicated expressions like:
Array<float,2> r = sum(A(i,k) * B(j,k), k);
How would this be done without expression templates? Would
A(i,k) * B(j,k) return a third-order tensor?
Another example from Boris's paper:
tensor tst1 = t2("ij")*t4("ijkl")*t4("klpq")*t2("pq");
This results in a scalar (zero-order tensor), since all the
indices are repeated.
If expression templates weren't used, can this still be
evaluated efficiently? (Hmm, can it be evaluated efficiently
even with expression templates??) The order of evaluation would have
a big effect on efficiency; depending on how you evaluate
it, you get 4th or 2nd order temporaries. Yuck, it is similar
to chains of matrix products.
--------------------- blitz-dev list --------------------------------
* To subscribe/unsubscribe: mail to majordomo@oonumerics.org, with
"subscribe blitz-dev" or "unsubscribe blitz-dev" in the body of the message
* Blitz++ web page: http://oonumerics.org/blitz/
This archive was generated by hypermail 2b29 : Wed Feb 20 2002 - 04:30:04 EST