|
|
Blitz Devel : |
From: jeremic_at_[hidden]
Date: 1998-07-17 16:04:19
Hello All,
On 14 Jul, Walter Landry wrote:
> 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.
> 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. 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.
> For example:
>
> tensor a = b("ijk") * c("klz") * d("z");
>
> a temorary will be created from c("klz") * d("z"), let's call it
> temp("kl") and then used in contraction with b("ijk") wich will create
> tensor a.
>
> Actually, I don't think a temp should be created. I think it should
> be inlined. My tensor package creates temps, though, because that was
> the only way I could think of doing it at the time. I'll discuss it
> more later.
That's the idea, use Expression Templates to get rid of temporaries (in
some cases).
>
>
> > If I wanted to contract an index internally, I would write
> >
> > Contraction_Index q('q);
> > S = T(q,q);
> >
>
> Or contract with Kronecker delta.
>
> But using a Kronecker delta seems as awkward as having to use a sum()
> function.
>
>
how about:
S = T.trace();
But that's matter of taste, one can put both or all three in there and
leave it to user to choose which to use.
> > K_new(ii,jj)<Lapse*(R(i,j) + Trace_K*K(i,j) - (2*K(i,k)&K_mix(k,j))
> > - 0.5*g(i,j)*(density - (g_up(k,l)|S(k,l))) - S(i,j));
> >
> In the above example I see only single tensorial multiplications, that
> is there is no tensor product with three tensors like in:
>
> tensor A = B("qwer") * c("erty") * D("tyui");
>
> or make it in Blitz notation:
>
> tensor A = B(q,,w,e,r) * c(e,r,t,y) * D(t,y,u,i);
>
> I'm not exactly sure what you're driving at, but I do have some
> examples that do this double multiplication
>
> 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.
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 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...
Best regards, Boris
Boris Jeremic Phone (315) 268-4435
Assistant Professor Fax (315) 268-7985
Department of Civil and Jeremic_at_[hidden]
Environmental Engineering http://www.clarkson.edu/~jeremic
Clarkson University
Potsdam, NY 13699-5710
--------------------- 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/