Hello All,
On 12 Jun, Todd L. Veldhuizen wrote:
> 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/
>
I have just updated the nDarray package on
http://civil.colorado.edu/nDarray/
so please download this latest version. It turns out that there was a
tiny memory leak (about 32 bytes every once in a while).
> 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..
>
I am particularly interested in computational mechanics, so my tensors
usually have 3 members in each dimension. So for example second order
tensor (stress tensor for example) would be a 3x3 matrix. People have
used tensors of 8-th order in theoretical developments (complex
hardening rules) but nobody has implemented anything of that sort
(think of using FORTRAN for something like that). I have set a variable
in nDarray that currently allows tensors up to 4 order mainly to make
it efficient to multiply them by parsing indices. But if somebody needs
higher ordet, it can be redefined.
>
> 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.
>
We can provide some basic tensors and leave it open for others to add
their own types. In nDaray all thensors are of the covariant type
although I have thought on implementing both upper and lower indeces
(covariant and contravariant). Only covariant tensors are used for
example in rectangular coordiante systems, while both covariant and
contraviariant are used in curvilinear coordinate systems.
I am further using nDarray tools for example in defining stresstensor
as a derived type from tensor (and straintensor and many others).
Input from other interested parties will be more than welcome!
>
> 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?
>
>
This turns out to be important if you have more than two tensors in
same command line. Consider:
tensor T = A("ij") * B("ijkl") * C("l");
here temporary resulting from A*B has to have it's own indices in order
to be properly parced and multiplied with C. If you use expression
templates you can unroll the whole line and then you probably don't
need to assign indices. But if there are more than three tensors in a
single command line, unrolling the expression might be tricky and
produce long code.
For nDarray tensor I usually set tensor indices to NULL between
operations just too prevent any suprises. For excample after previos
command line I might say:
T.null_indices();
and be sure that in next use of tensor T it will start with fresh
indices. This is also done at the end of all operators (operator*,
operator/ . . . ) and I am just doing again to be sure about having
fresh indices every time I need them.
> 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".
>
With nDarray you can contract more than one index on the same command
line. This is usually needed in computational mechanics, for example
tensor stress = E("ijkl") * epsilon("kl");
will return proper tensor for stress from contracting kl indices from
stiffness tensor (fourth order) E and strain tensor epsilon.
>
> 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.
>
IN this case there are couple of additional concerns. This line of code
actually is not portable. Following is from the new math_tst.cc on the
above mentioned web site:
// Some relevant notes from the Dec. 96 Draft (as pointed out by KAI folks):
//
// "A binary operator shall be implemented either by a non-static member
// function (_class.mfct_) with one parameter or by a non-member function
// with two parameters. Thus, for any binary operator @, x@y can be
// interpreted as either x.operator@(y) or operator@(x,y)." [13.5.2]
//
// "The order of evaluation of arguments [of a function call] is
// unspecified." [5.2.2 paragraph 9]
//
//
// SO THAT IS WHY I NEED THIS temp
// It is compiler dependent (bad) which operator will be called up first:
// operator() or operator* , so define tensor temp. Needed only if
// the same tensor appears more than once in the same command line.
// for example this is compiler dependent (with some compilers it works)
// tst1 = t2("ij") * t4("ijkl") * t4("klpq") * t2("pq");
// so it's better to use tensor temp:
tensor tst1;
tensor temp = t2("ij") * t4("ijkl");
tst1 = temp("kl") * t4("klpq") * t2("pq");
tst1.print("tst1","\ntensor tst1 = t2(\"ij\")*t4(\"ijkl\")*t4(\"klpq\")*t2(\"pq\");");
It turns out that the line
tst1 = t2("ij") * t4("ijkl") * t4("klpq") * t2("pq");
works with compilers that evaluate operator* first and doesn't work
with others. Putting brackets doesn't help with the order of evaluation
of operator* or operator(). It only helps if you want to control
efficiency of creation of temporaries.
So much for now. Comments are more than welcome!
Best reagards, Boris
Boris Jeremic
Assistant Professor
Rowley Labs 236
Department of Civil and Environmental Engineering
Clarkson University
Potsdam, NY 13699-5710
Phone (315) 268-4435
Fax (315) 268-7985
Jeremic@Polaris.Clarkson.edu
http://www.clarkson.edu/~jeremic
--------------------- 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