|
|
Blitz Devel : |
From: Walter Landry (landry_at_[hidden])
Date: 1998-07-20 17:25:10
> I think that indices should ALWAYS be required. I don't consider it a
> hassle to put them in, and it makes the code much more readable. I
> think that forcing people to write good code is a good thing, as long
> as you don't lose any power in the process. Also, looking at your
you will be invoking an overloaded operator() for no reason.
With Expression Templates, I don't think you have to invoke anything
extra at all. The expressions are going to get unrolled anyway.
> example, you have the statement
>
> 2.*K("ik")*K_mix("kj")
>
> whereas I have
>
> 2*K(i,k)&K_mix(k,j)
>
> In this example, K(i,k) is a symmetric rank 2 tensor, and K_mix(k,j)
> is not symmetric. I use the "&" symbol because I want to produce a
> rank 2 symmetric matrix out of it. This only requires computation of
> 6 elements. I have a feeling that your routines compute all 9. There
> has to be a way to tell the compiler that you only want 6 elements, or
> efficiency will suffer *greatly*. I don't know how to do this except
> by using a different symbol for the contraction (which is what I do).
>
It will not really suffer greatly. Difference is 33% but you are right
it might be worth looking into making special types.
The difference is greater when you start looking at antisymmetric
tensors. For example, the Riemann (a common tensor in General
Relativity) has four indices and, in three dimensions, only six
independent components. There the difference is more like a factor of
ten. In any case, I don't think we should be throwing away
efficiencies of 33%. Decisions like that make people not use the
tensor package.
As to the above expression, how about:
2.*(K("ik")*K_mix("kj")).simmetrize();
or symmetrize(2.*K("ik")*K_mix("kj"));
This of course works only for 2 tensor. For 4 tensor you have to
specify which symmetry you want, major, minor or some combination.
This may turn into a matter of taste, but I must say that I prefer
using funny symbols to mean add or contract as opposed to writing
things like symmetrize(). It is more compact. Neither are good
solutions. We might be able to write it so that it looks at what the
eventual result is. For example, look at
tensor K_new = Lapse*R("ij") + Trace_K*K("ij") - (2.*K("ik")*K_mix("kj") -
0.5 * g("ij") * (density - g_up("kl")*S("kl")) - S("ij");
K_new is a symmetric rank 2 tensor. We might be able to get it to
parse the expression so that it will "know" that when
(2.*K("ik")*K_mix("kj")
is calculated, it will be going into a symmetric tensor. Then it will
only calulate the six components. We might even get this for free,
since when the expression is unrolled, it will only make six
expressions, not nine. This approach has the advantage of being most
like what people do when they look at an equation (they know that the
result is symmetric).
However, then the computer doesn't do any type checking. People do.
All rank 2 tensors (dense, symmetric, antisymmetric) are freely
interchangeable. That may not be such a bad thing, since it is not
unreasonable to expect people to debug their formulas before entering
them into the computer.
Also, I have a nagging suspicion that there will be cases where it
won't work (it will compute more than it should). In my own
complicated code, I haven't found any examples, but I could make one
up:
Tensor0 W = K(i,k)*K_mix(k,j)*g(i,j);
Using the unrolling/temp rules I proposed in an earlier email, this
would make a temp from K(i,k)*K_mix(k,j), that would then be
contracted with g(i,j) (it evaluates left to right). If we don't have
a special contraction operator, it becomes 9+9=18 multiplies. A
little worse than the 6+9=15 with a special contraction operator (I
really don't know if I am measuring the most relevant thing here).
I guess the upshot is that it may be worth it to have a special
operator (whether it is another overloaded operator, or a symmetrize()
function) to cover cases like this, but it won't be used that much.
If it is rare, it doesn't matter as much if it is ugly. I'm not so
sure that it is rare, though. It's not in my code, but that doesn't
mean all that much.
> lX(ii,jj) < (g_up(i,k)%dX(j,k) || g_up(j,k)%dX(i,k))
> - (2/3.0)*dX(q,q)*g_up(i,j)
> + (g_up(i,k)%(christof(j,k,l)%X(l)) || g_up(j,k)%(christof(i,k,l)%X(l)));
>
> The "||" symbol adds two non_symmetric rank 2 tensors to make a
> symmetric rank 2 tensor. Once again, it only does 6 instead of 9
> adds.
>
There is a problem here since C++ standard does not say what the order
of execution of functions should be. The above line might not be
portable across different compilers. I am not sure about the details of
your code, but try compiling and runing with different compilers.
Are you saying that C++ is getting rid of operator precedence? That
is all that I use to uniquely define how everything above is
evaluated. That is why I have parentheses in varying places.
I am preaty busy with some stuff I need to finish in next 10 or so
days, but I'd like to continue this discussion. August might be good
month to actually start some work on tensors (combine existing
developments and look for improvements).
I am going on vacation during the second and third week of August, and
I'm not sure in general how much time I can spend on it. It depends
on how my thesis is coming.
I would appreciate if others interested in tensors would raise their
voice as to what they want to see in a tensor package or what is their
application...
As would I.
Sincerely,
Walter Landry
landry_at_[hidden]
--------------------- blitz-dev list --------------------------------
* To subscribe/unsubscribe: mail to majordomo_at_[hidden], with
"subscribe blitz-dev" or "unsubscribe blitz-dev" in the body of the message
* Blitz++ web page: http://oonumerics.org/blitz/