## digitalmars.D - Thoghts on: opSlice, array expressions & vectorization

- Norbert Nemec <Norbert.Nemec gmx.de> Jul 01 2004
- Patrick Down <Patrick_member pathlink.com> Jul 02 2004
- Daniel Horn <hellcatv hotmail.com> Jul 02 2004
- Norbert Nemec <Norbert.Nemec gmx.de> Jul 02 2004
- Regan Heath <regan netwin.co.nz> Jul 02 2004
- Patrick Down <pat codemoon.com> Jul 02 2004
- "Ivan Senji" <ivan.senji public.srce.hr> Jul 02 2004
- Norbert Nemec <Norbert.Nemec gmx.de> Jul 02 2004
- Stephen Waits <steve waits.net> Jul 02 2004
- Norbert Nemec <Norbert.Nemec gmx.de> Jul 03 2004
- =?iso-8859-1?q?Knud_S=F8rensen?= <knud NetRunner.all-technology.com> Jul 03 2004
- Norbert Nemec <Norbert.Nemec gmx.de> Jul 03 2004
- "Dan Williams" <dnews ithium.NOSPAM.net> Jul 03 2004
- =?iso-8859-1?q?Knud_S=F8rensen?= <knud NetRunner.all-technology.com> Jul 03 2004
- Norbert Nemec <Norbert.Nemec gmx.de> Jul 03 2004
- Stephen Waits <steve waits.net> Jul 03 2004
- Norbert Nemec <Norbert.Nemec gmx.de> Jul 04 2004
- Stephen Waits <steve waits.net> Jul 05 2004
- Norbert Nemec <Norbert.Nemec gmx.de> Jul 06 2004
- =?iso-8859-1?q?Knud_S=F8rensen?= <knud NetRunner.all-technology.com> Jul 02 2004
- Norbert Nemec <Norbert.Nemec gmx.de> Jul 03 2004
- "Walter" <newshound digitalmars.com> Jul 02 2004
- Norbert Nemec <Norbert.Nemec gmx.de> Jul 03 2004
- "Ivan Senji" <ivan.senji public.srce.hr> Jul 03 2004
- Norbert Nemec <Norbert.Nemec gmx.de> Jul 03 2004
- "Ivan Senji" <ivan.senji public.srce.hr> Jul 03 2004
- "Dan Williams" <dnews ithium.NOSPAM.net> Jul 03 2004
- Stephen Waits <steve waits.net> Jul 05 2004
- Norbert Nemec <Norbert.Nemec gmx.de> Jul 05 2004
- Andy Friesen <andy ikagames.com> Jul 03 2004
- Norbert Nemec <Norbert.Nemec gmx.de> Jul 03 2004

Hi there, just a few words about this nights thoughts. This is probably post-1.0, but it should already be considered now in order to avoid steps in the wrong direction. For high-performance numerics, vectorization of array operations is absolutely crucial. Currently, I know of several approaches to this problem in different languages: * HPF (High-Performance-Fortran) has powerful vector expression capabilities built in. AFAIK, these are always limited to one routine, i.e. you cannot spread vectorization across routines. This makes it mostly impossible to extend the existing set of vector operations. * Some very promising languages (e.g. Single-Assignment-C: http://www.sac-home.org) take the step to purely functional languages, or even purely declarative languages, where the compiler can optimize much more aggressively than in imperative languages. * C++ has expression templates which allow to resolve the vectorization of expression to be done via meta-programming. This allows utmost flexibility, but is very clumsy to use, and what is even more problematic: it can never be tuned as closely to the processor as a good optimizing compiler can. For D, I already have some ideas how these problems could be overcome, finding a way in the middle. -> Full support by the language -> Simple to use -> Simple to extend -> compatible with the (mostly) imperative environment of D Core idea is, to drop the current concept of array expressions in favor of a more flexible basic concept of vector expressions, using index notation. (I'm not sure about the syntax yet, so I cannot give conclusive details) The functionality of the current (specified but not implemented) array expressions could then simply be implemented in the library without performance penalty. The other pillar on which the whole idea rests is slice assignment. Thinking about the matter for some time, I think, it might not be a good idea to make slice assignments overloadable at the moment. Having opIndex and opIndexAssign is a nice, unproblematic thing: as soon as you index an array, you break vectorization anyway, so any fully user-defined array may perform just as good as native arrays. Handling slices though, is tightly coupled with vectorizing code. Once, we know how vector expressions really work, we may find some way to offer a overloadable opSlice and maybe even opSliceAssign. Until then, we should better leave slicing to the native arrays as we have them and not do any rash decisions trying to force multidimensional slicing into the language without having multidimensional arrays. Even the existing opSlice should probably be dropped for now. It certainly is easier to reintroduce later when we realize that it fits in than to remove it if it collides with the rest of the vectorization/slicing concept. OK, so far for now. Actually, this post mostly means that the matters are postponed further, but it also means that I'm still actively thinking about them. Ciao, Nobbi

Jul 01 2004

In article <cc2un6$11gc$1 digitaldaemon.com>, Norbert Nemec says...Hi there, just a few words about this nights thoughts. This is probably post-1.0, but it should already be considered now in order to avoid steps in the wrong direction.

Let me throw out a thought here. How would you feel about an explicit vectorblock keyword? Inside a vectorblock the compiler has a restricted/exteneded subset of the language that can be aggressively optimized for vector operatons. What I worry about is that we may hurt the general expressiveness of the language if we put too many restrictions just for vector operations. I am already annoyed that I can't do things like a[4..12] = a[3..11]. But inside a vectorblock you can make assumptions about array aliasing and other things. The compiler would insert code in debug mode to check these assuptions at the begining of the block.

Jul 02 2004

It's beginning to sound like my project (in C though) called brook http://brook.sf.net/ still has a lot of work to become useful on cpu's or on clusters...but similar idea...:-) Patrick Down wrote:In article <cc2un6$11gc$1 digitaldaemon.com>, Norbert Nemec says...Hi there, just a few words about this nights thoughts. This is probably post-1.0, but it should already be considered now in order to avoid steps in the wrong direction.

Let me throw out a thought here. How would you feel about an explicit vectorblock keyword? Inside a vectorblock the compiler has a restricted/exteneded subset of the language that can be aggressively optimized for vector operatons. What I worry about is that we may hurt the general expressiveness of the language if we put too many restrictions just for vector operations. I am already annoyed that I can't do things like a[4..12] = a[3..11]. But inside a vectorblock you can make assumptions about array aliasing and other things. The compiler would insert code in debug mode to check these assuptions at the begining of the block.

Jul 02 2004

Patrick Down wrote:In article <cc2un6$11gc$1 digitaldaemon.com>, Norbert Nemec says...Hi there, just a few words about this nights thoughts. This is probably post-1.0, but it should already be considered now in order to avoid steps in the wrong direction.

Let me throw out a thought here. How would you feel about an explicit vectorblock keyword? Inside a vectorblock the compiler has a restricted/exteneded subset of the language that can be aggressively optimized for vector operatons. What I worry about is that we may hurt the general expressiveness of the language if we put too many restrictions just for vector operations. I am already annoyed that I can't do things like a[4..12] = a[3..11]. But inside a vectorblock you can make assumptions about array aliasing and other things. The compiler would insert code in debug mode to check these assuptions at the begining of the block.

First: of course, D is a general purpose language. Trading general usefulness for qualities in some special area is, of course, not an option. Second: I don't really like the idea of splitting the language into different "modes". Once you start splitting, the parts will begin to drift apart and in the end you have a mess of different languages. Conclusion: Of course, vector operations should not restrict the language, but they should still be built into the language with the optimum combination of performance and practicability. Now your example: You are stating a vector-expression. One possible definition for a vector expression is that the order of execution of the individual parts is not defined. The effect of this is, that your statement is not an error (the compiler cannot detect it), but the outcome is not defined. Some possible ideas that you might have had to solve the problem would be: * either define a fixed order of execution in some arbitrary way by the language specs - this would be just cause confusion without solving the problem. (Suddenly a[3..11] = a[4..12] would be correct, but a[4..12] = a[3..11] would be incorrect!) * ask for enough intelligence in the compiler to choose the right order - I'm pretty sure this is impossible since aliasing cannot be detected at compile time. Your example seems simple enough, but the general case would be a mess to sort out. * introduce temporary copies wherever aliasing might be possible - this would mean *lots* of overhead that does not only hurt high-performance computing but everyone. There might be other approaches, but I doubt there is a simple solution. Of course, all of this has to be investigated more closely, but I don't think anyone would come up with a good solution that allows statements like yours. Ciao, Nobbi

Jul 02 2004

On Fri, 2 Jul 2004 16:24:06 +0000 (UTC), Patrick Down <Patrick_member pathlink.com> wrote:

Let me throw out a thought here. How would you feel about an explicit vectorblock keyword? Inside a vectorblock the compiler has a restricted/exteneded subset of the language that can be aggressively optimized for vector operatons. What I worry about is that we may hurt the general expressiveness of the language if we put too many restrictions just for vector operations. I am already annoyed that I can't do things like a[4..12] = a[3..11]. But inside a vectorblock you can make assumptions about array aliasing and other things. The compiler would insert code in debug mode to check these assuptions at the begining of the block.

What do you want a[4..12] = a[3..11] to do exactly? It looks to me like you want to basically do a memmove(&a[4],&a[3],typeof(a[0]).sizeof * 8); or is it supposed to do something else? Regan -- Using M2, Opera's revolutionary e-mail client: http://www.opera.com/m2/

Jul 02 2004

Regan Heath <regan netwin.co.nz> wrote in news:opsajaypgc5a2sq9 digitalmars.com:What do you want a[4..12] = a[3..11] to do exactly? It looks to me like you want to basically do a memmove(&a[4],&a[3],typeof(a[0]).sizeof * 8);

Yes, exactly what b[4..12] = a[3..11] does right now in D. a[4..12] = a[3..11] won't work because the memory areas overlap.

Jul 02 2004

I am pretty interested in the future improvements of D's arrays, so can you please explain (in a couple of words) what is vectorization (or point me to some place where i could find out more). As for array operations: The only thing that is currently implemented is int[] x= new int[100]; x[]=3; and this maybe even isn't an array operation? But the way i see this feature is like just another way of writing a loop on the elements of the array. For example z[] = x[]+b[]; i see as: for(int i=0;iyz.length;i++) z[i]=x[i]+b[i]; So i think this feature could be more general by maybe introducing a new syntax to say explicitly "for every element of the array" for example: int[] x,y,z; z[#] = x[#]+y[#]; And this would mean mathematically speaking z[i]=x[i]+y[i]; We could even do something like: int func(int x); func(z[#]); //call func with every z[i], z[i] is int and func takes int so everything is ok :) Access to the index would be even greater. z[#(int index1)] = 5*index1-3; or even int[10][10] multidim; multidim[#(int i)][#(int j)] = i*j+1; //some matrix initialization or multidim[#(int i)] = z; //meaning "for every posible first index assign to multidim[i] = z" Quite possible many of the things i said are even impossible to implement or maybe useless, but i do beleive array operations could be made very powerfull and usefull feature (maybe not even in this direction but in some other) Norbert, tell me if i am totally crazy :)

Jul 02 2004

Ivan Senji wrote:I am pretty interested in the future improvements of D's arrays, so can you please explain (in a couple of words) what is vectorization (or point me to some place where i could find out more)

Norbert, tell me if i am totally crazy :)

No, your thoughts go pretty much the right direction. But as you can see yourself, it really is difficult to find a syntax that is powerful, clean, simple and extensible. As always in programming language design, you want to have the simple things simple and the complex ones possible. I have not yet come up with a complete concept for the syntax, but there are some ideas for most pieces already. One of the core ideas is that you would need two different ways for writing vector expressions. Index-free like in a[] = b[1..3] + c[2..4]; for the simple cases. Here you would need to slice all the operands into the right shape and then add then. An alternative way to specify the same thing would be something like: a[] = [i=1..3](b[i] + c[i+1]); This is, of course, much more flexible, also allowing complicated computation inside the vector statement. In the end, you could of course combine both ways wildly. One important detail is, that vectorization should be possible across function calls. That is, if you write and use something like int[] myfunc(int[] a,int[] b) { return a + b; } the language should be defined in a way that allows the compiler to inline this without breaking vectorization. (i.e. without creating intermediate temporary arrays)

Jul 02 2004

Norbert Nemec wrote:OK, so far for now. Actually, this post mostly means that the matters are postponed further, but it also means that I'm still actively thinking about them.

Great post Norbert.. I'd _LOVE_ to see something like this brought into the language. It opens up an avenue for D in the scientific and real-time-3D area which, IMO, does not currently exist. My only suggestion is that these vector expressions be able to expand to handle non-trivial cases. For example, think of calculating a vector cross product, or matrix multiplication, quaternion multiplication, etc. What I desire is some of the performance of Boost, without the HORRID ugliness. In C++ this is accomplished through template expressions, and template meta-programming, all which sounds fun - but in practice it is just awful to work with. This is all an effort to escape from TEMPORARY OBJECT HELL.. it's an especially hellish hell. Optimization opportunities that compilers like Codeplay's take advantage of would elevate this language WAY up there IMO.. I believe this language CAN be built to make this easier on compilers, WITHOUT compromising the general purpose goals. Ideally, something like: Vector3 a, b, c; c = 2.0 * a + b; Should reduce and be inlined to something like: <vectorizable> c[0] = 2.0 * a[0] + b[0]; c[1] = 2.0 * a[1] + b[1]; c[2] = 2.0 * a[2] + b[2]; </vectorizable> Or, a bigger, more complicated example: Vector3 a, b, c; Matrix33 m; c = m * a.Cross(b); becomes: <vectorizable> c[0] = a[0]*(b[1]*m[0,2] - b[2]*m[0,1]) + a[1]*(b[2]*m[0,0] - b[0]*m[0,2]) + a[2]*(b[0]*m[0,1] - b[1]*m[0,0]); c[1] = a[0]*(b[1]*m[1,2] - b[2]*m[1,1]) + a[1]*(b[2]*m[1,0] - b[0]*m[1,2]) + a[2]*(b[0]*m[1,1] - b[1]*m[1,0]); c[2] = a[0]*(b[1]*m[2,2] - b[2]*m[2,1]) + a[1]*(b[2]*m[2,0] - b[0]*m[2,2]) + a[2]*(b[0]*m[2,1] - b[1]*m[2,0]); </vectorizable> --Steve PS - Your other post (below) about aliasing is also very important. In game development, aliasing is one of our worst enemies.

Jul 02 2004

Stephen Waits wrote:My only suggestion is that these vector expressions be able to expand to handle non-trivial cases. For example, think of calculating a vector cross product, or matrix multiplication, quaternion multiplication, etc.

I don't think, that would be a good idea: Plain vector operations cover *elementwise* operations. With a good concept, it should be possible to stay flexible and extensible enough, if this is put into the language. The operations you are suggesting go far beyond that in complexity. Of course, you could great performance having matrix multiplications etc. built into the language and optimized by the compiler, but that would be far beyond what an extensible general purpose language should do. The better way to go here would be to go for expression template programming which should just be handled in a cleaner way than in C++. There, you can implement your linear/quaternion or whatever algebra and hand-tune it for your needs. Vector expressions would certainly help you some way, since they already eliminate much of the problems involving temporary objects, but beyond that meta-programming would be the way to go. Ciao, Nobbi PS: Since you are talking about quaternions - have you ever heard about Geometric Algebra (as introduced by David Hestenes)? For handling geometry in mathematics and computation, this is even more powerful than quaternions. Takes a while to get used to that kind of mathematical language, but once you've understood the key concepts, geometry suddenly becomes so much clearer! (And code possibly more efficient as well - there is some promising research going on in that direction.)

Jul 03 2004

On Sat, 03 Jul 2004 09:51:11 +0200, Norbert Nemec wrote:Stephen Waits wrote:My only suggestion is that these vector expressions be able to expand to handle non-trivial cases. For example, think of calculating a vector cross product, or matrix multiplication, quaternion multiplication, etc.

I don't think, that would be a good idea: Plain vector operations cover *elementwise* operations. With a good concept, it should be possible to stay flexible and extensible enough, if this is put into the language.

Einstein notation would give you all this and much more. it would also give you powerful vectorization and avoid temporaries. agian I invite you to take a look at http://dlanguage.netunify.com/56 and tell me what you think. Knud

Jul 03 2004

Knud Sørensen wrote:On Sat, 03 Jul 2004 09:51:11 +0200, Norbert Nemec wrote:Stephen Waits wrote:My only suggestion is that these vector expressions be able to expand to handle non-trivial cases. For example, think of calculating a vector cross product, or matrix multiplication, quaternion multiplication, etc.

I don't think, that would be a good idea: Plain vector operations cover *elementwise* operations. With a good concept, it should be possible to stay flexible and extensible enough, if this is put into the language.

Einstein notation would give you all this and much more. it would also give you powerful vectorization and avoid temporaries. agian I invite you to take a look at http://dlanguage.netunify.com/56 and tell me what you think.

See my other message about it. Just one point for clarification: Something like: a[i,j] = sum(k){b[i,k]*c[k,j]}; // Pseudo-syntax! is called *index-notation* - this is exactly the core of my concept for vector operations. (I have no idea yet how to do aggregators like "sum", but they will definitely be needed) It was used in mathematics long before Einstein was born. The term "Einstein-Notation" on the other hand only refers to dropping the "sum" in the above expression, saying this sum is implicitly to be understood if one index appears twice. Einstein introduced that convention because he was lazy and because it really saves on writing in relativistic calculations. In contexts outside of relativity, geometry and linear algebra, this convention does not make much sense, though and actually even gets in the way. Therefore: index-notation makes perfect sense in D, Einstein-notation (i.e. summing over indices implicitly) would just be extremely confusing and counterproductive in many areas.

Jul 03 2004

"Norbert Nemec" <Norbert.Nemec gmx.de> wrote in message news:cc5t00$2m4h$1 digitaldaemon.com...Knud Sørensen wrote:On Sat, 03 Jul 2004 09:51:11 +0200, Norbert Nemec wrote:Stephen Waits wrote:My only suggestion is that these vector expressions be able to expand

handle non-trivial cases. For example, think of calculating a vector cross product, or matrix multiplication, quaternion multiplication,

I don't think, that would be a good idea: Plain vector operations cover *elementwise* operations. With a good concept, it should be possible to stay flexible and extensible enough, if this is put into the language.

Einstein notation would give you all this and much more. it would also give you powerful vectorization and avoid temporaries. agian I invite you to take a look at http://dlanguage.netunify.com/56 and tell me what you think.

See my other message about it. Just one point for clarification: Something like: a[i,j] = sum(k){b[i,k]*c[k,j]}; // Pseudo-syntax! is called *index-notation* - this is exactly the core of my concept for vector operations. (I have no idea yet how to do aggregators like "sum", but they will definitely be needed) It was used in mathematics long before Einstein was born. The term "Einstein-Notation" on the other hand only refers to dropping the "sum" in the above expression, saying this sum is implicitly to be understood if one index appears twice. Einstein introduced that convention because he was lazy and because it really saves on writing in relativistic calculations. In contexts outside of relativity, geometry and linear algebra, this convention does not make much sense, though and actually

gets in the way. Therefore: index-notation makes perfect sense in D, Einstein-notation

summing over indices implicitly) would just be extremely confusing and counterproductive in many areas.

I'm with Norman on this... in a programming langauge there has to be as much expliciticity (is that a word??!?) and as little ambiguity as possible. I feel that Einstein notation would quickly get confusing because of the lack of explicit operation identifiers - methodology by context is a nice idea, but impractical for this I believe. Also, wouldn't it have the side-effect of being more work for the compiler? Whereas index notation would not have that disadvantage.

Jul 03 2004

On Sat, 03 Jul 2004 11:07:50 +0200, Norbert Nemec wrote:Knud Sørensen wrote:On Sat, 03 Jul 2004 09:51:11 +0200, Norbert Nemec wrote:Stephen Waits wrote:

I don't think, that would be a good idea: Plain vector operations cover *elementwise* operations. With a good concept, it should be possible to stay flexible and extensible enough, if this is put into the language.

Einstein notation would give you all this and much more. it would also give you powerful vectorization and avoid temporaries. agian I invite you to take a look at http://dlanguage.netunify.com/56 and tell me what you think.

See my other message about it. Just one point for clarification:

Yes, we did post replys to each other at the same time.Something like: a[i,j] = sum(k){b[i,k]*c[k,j]}; // Pseudo-syntax! is called *index-notation* - this is exactly the core of my concept for vector operations. (I have no idea yet how to do aggregators like "sum", but they will definitely be needed) It was used in mathematics long before Einstein was born.

Yes, but how would you implement this sum aggregator on an expression like # B[i,j] = A[k,i]*C[k,l]*A[l,j]; // B=A^T*C*A; without devectorization. What Einstein-Notation says that we can drop the sum because the order of the sums don't matter. It is in fact a vectorization of multi summations.The term "Einstein-Notation" on the other hand only refers to dropping the "sum" in the above expression, saying this sum is implicitly to be understood if one index appears twice. Einstein introduced that convention because he was lazy and because it really saves on writing in relativistic calculations. In contexts outside of relativity, geometry and linear algebra, this convention does not make much sense, though and actually even gets in the way. Therefore: index-notation makes perfect sense in D, Einstein-notation (i.e. summing over indices implicitly) would just be extremely confusing and counterproductive in many areas.

Can you give an example of an vectorization which get extremely confusing and counterproductive in Einstein notation. What about making a vectorization notation shootout. It could be a wiki page where one could post vectorization problems, an we could post solutions in different notations. Then when enough different problems/solutions where posted it would be easier to compare the notations. Knud

Jul 03 2004

Knud Sørensen wrote:Something like: a[i,j] = sum(k){b[i,k]*c[k,j]}; // Pseudo-syntax! is called *index-notation* - this is exactly the core of my concept for vector operations. (I have no idea yet how to do aggregators like "sum", but they will definitely be needed) It was used in mathematics long before Einstein was born.

Yes, but how would you implement this sum aggregator on an expression like # B[i,j] = A[k,i]*C[k,l]*A[l,j]; // B=A^T*C*A; without devectorization.

Simple: B[i,j] = sum(k){sum(l){A[k,i]*C[k,l]*A[l,j]}}; or, equivalently: B[i,j] = sum(k){A[k,i]*sum(l){C[k,l]*A[l,j]}}; B[i,j] = sum(l){sum(k){A[k,i]*C[k,l]*A[l,j]}}; B[i,j] = sum(l){sum(k){A[k,i]*C[k,l]}*A[l,j]}; B[i,j] = sum(l){A[l,j]*sum(k){A[k,i]*C[k,l]}}; B[i,j] = sum(l){sum(k){A[k,i]*C[k,l]*A[l,j]}};What Einstein-Notation says that we can drop the sum because the order of the sums don't matter. It is in fact a vectorization of multi summations.

This is wrong. It is a general fact that you can switch the order of the sums and that you can pull out factors that do not depend on the summation index. Einstein notation is only possible because of this fact. But for every Einstein expression, you can simply put in the summation signs in some appropriate position. Really: Einstein notation is *always* just a shorthand for an explicit summation.Can you give an example of an vectorization which get extremely confusing and counterproductive in Einstein notation.

D[i] = sum(l){A[l,i]*B[l,i]*C[l,i]}; D[i,j] = A[i,j]*B[i,j]; // no summation here - just elementwise // multiplication As said before: these do not make sense for matrices (in the linear algebra meaning), but for general arrays, they may be perfectly sensible.What about making a vectorization notation shootout. It could be a wiki page where one could post vectorization problems, an we could post solutions in different notations.

I don't think that would lead very far. What would be needed is a complete concept for arrays and vectorization. Notation only follows from this concept. If there are several competing concepts we can start discussing. So far, I don't even have a full concept myself. This really is not the time for discussing notation.

Jul 03 2004

On Sat, 03 Jul 2004 18:45:10 +0200, Norbert Nemec wrote:Knud Sørensen wrote:Something like: a[i,j] = sum(k){b[i,k]*c[k,j]}; // Pseudo-syntax! is called *index-notation* - this is exactly the core of my concept for vector operations. (I have no idea yet how to do aggregators like "sum", but they will definitely be needed) It was used in mathematics long before Einstein was born.

Yes, but how would you implement this sum aggregator on an expression like # B[i,j] = A[k,i]*C[k,l]*A[l,j]; // B=A^T*C*A; without devectorization.

Simple: B[i,j] = sum(k){sum(l){A[k,i]*C[k,l]*A[l,j]}}; or, equivalently: B[i,j] = sum(k){A[k,i]*sum(l){C[k,l]*A[l,j]}}; B[i,j] = sum(l){sum(k){A[k,i]*C[k,l]*A[l,j]}}; B[i,j] = sum(l){sum(k){A[k,i]*C[k,l]}*A[l,j]}; B[i,j] = sum(l){A[l,j]*sum(k){A[k,i]*C[k,l]}}; B[i,j] = sum(l){sum(k){A[k,i]*C[k,l]*A[l,j]}};

My question where not about notation but about implimention. My fear is that the implimention would do the devectorization.What Einstein-Notation says that we can drop the sum because the order of the sums don't matter. It is in fact a vectorization of multi summations.

This is wrong. It is a general fact that you can switch the order of the sums and that you can pull out factors that do not depend on the summation index. Einstein notation is only possible because of this fact. But for every Einstein expression, you can simply put in the summation signs in some appropriate position. Really: Einstein notation is *always* just a shorthand for an explicit summation.Can you give an example of an vectorization which get extremely confusing and counterproductive in Einstein notation.

D[i] = sum(l){A[l,i]*B[l,i]*C[l,i]}; D[i,j] = A[i,j]*B[i,j]; // no summation here - just elementwise // multiplication

With Einstein notation you would write # D[i] = A[l,i]*B[l,i]*C[l,i]; Index l is a free index in a multiplication, so there is a summation. # D[i,j] = A[i,j]*B[i,j]; There is no free index, so the is no summation. A free index is a index which is only defined on the right hand side.As said before: these do not make sense for matrices (in the linear algebra meaning), but for general arrays, they may be perfectly sensible.What about making a vectorization notation shootout. It could be a wiki page where one could post vectorization problems, an we could post solutions in different notations.

I don't think that would lead very far. What would be needed is a complete concept for arrays and vectorization. Notation only follows from this concept. If there are several competing concepts we can start discussing. So far, I don't even have a full concept myself. This really is not the time for discussing notation.

Yes, but how do you know that the concept is complete. If you don't have a set of problems to compare it with ?? It is also better to have the problems first or else one tend to find problems which fit nice into the concept, instead of a concept that fits the actual problems.

Jul 03 2004

Knud Sørensen wrote:My question where not about notation but about implemention. My fear is that the implemention would do the devectorization.

That depends on how much effort your put into the optimization. Vector expressions are just a means to *allow* the compiler to optimize. Of course, once this is given, a lot of effort and research has to be put into the implementation to really squeeze out that performance. A simple (and correct) implementation would be to just turn the vector expression into a loop. There would already be some gain compared to traditional loop-programming, since this vectorization can be threaded through function calls, avoiding unnecessary temporary copies of the data. More sophisticated implementations could go a lot further, sorting the commands for optimal cache/pipeline usage, or maybe even distributed execution. A language never is fast or slow by itself, but different languages make it easier or harder for the compiler to optimize the code.What Einstein-Notation says that we can drop the sum because the order of the sums don't matter. It is in fact a vectorization of multi summations.

This is wrong. It is a general fact that you can switch the order of the sums and that you can pull out factors that do not depend on the summation index. Einstein notation is only possible because of this fact. But for every Einstein expression, you can simply put in the summation signs in some appropriate position. Really: Einstein notation is *always* just a shorthand for an explicit summation.Can you give an example of an vectorization which get extremely confusing and counterproductive in Einstein notation.

D[i] = sum(l){A[l,i]*B[l,i]*C[l,i]}; D[i,j] = A[i,j]*B[i,j]; // no summation here - just elementwise // multiplication

With Einstein notation you would write # D[i] = A[l,i]*B[l,i]*C[l,i]; Index l is a free index in a multiplication, so there is a summation.

But here you will still have to say that l is a index at all. As it stands, it is just an undeclared variable name. So you can either construct some syntax like index l; D[i] = A[l,i]*B[l,i]*C[l,i]; or stick with simply do some syntax like: D[i] = sum(l){A[l,i]*B[l,i]*C[l,i]}; Which is much more descriptive. B.t.w: one other example where you will run into trouble: D[i] = sin(A[l,i]*B[l,i]); Now: should that implicit summation happen inside or outside the sine? How about: D[i] = sin(A[l,i])*sin(B[l,i]); And, what about other accumulators? There should be at least an accumulative "product", as well as "and" and "or" as primitives. (Maybe others as commodity: "mean" or "stdev" etc.)A free index is a index which is only defined on the right hand side.

I understand pretty well, but I still don't think that would be a good idea in a programming language.Yes, but how do you know that the concept is complete. If you don't have a set of problems to compare it with ??

I'm not talking of "complete" in some mathematical sense. In this context I would call a "complete" concept one that can be written down in a form that is fit for the language specs. The behavior has to be clearly defined not only for examples but for the general case.It is also better to have the problems first or else one tend to find problems which fit nice into the concept, instead of a concept that fits the actual problems.

Sometimes, it might yet be the other way around. A good programming language does not only solve problems but it inspires you to solve problems you would never have thought of before. Anyway: in this case, my ideas came from thinking about a very real problem. Vectorization is the key to performance on modern hardware. C++ specialists went through the pains of writing the boost++ library to solve this problem. I think it is worth thinking about ways of doing better in D. Furthermore: I've been using Matlab for my numerics work lately and I've come to see how much a powerful array mechanism helps in programming. Some of my recent ideas were certainly inspired by Matlab's strengths, others by its weaknesses.

Jul 03 2004

Norbert Nemec wrote:Stephen Waits wrote:

I don't think, that would be a good idea: Plain vector operations cover *elementwise* operations. With a good concept, it should be possible to stay flexible and extensible enough, if this is put into the language. The operations you are suggesting go far beyond that in complexity. Of

Sorry, I wasn't suggesting that these operations be built into the language; only that the language allow enough flexibility so that we may code such optimizations without sacrificing so much performance. I shall research the geometric algebra you mention. We generally just use quaternions for character animation, and possibly camera animation. --Steve

Jul 03 2004

Stephen Waits wrote:Sorry, I wasn't suggesting that these operations be built into the language; only that the language allow enough flexibility so that we may code such optimizations without sacrificing so much performance.

OK, sorry. If that is so, we are completely on line there.I shall research the geometric algebra you mention. We generally just use quaternions for character animation, and possibly camera animation.

Well, quaternions are a convenient way to express rotations in 3D-space, just like complex numbers are very convenient for 2D rotations. Both are more efficient and simpler than the usual matrix-algebra. Geometric algebra just takes this to a more general level working for spaces of arbitrary dimension (also relativistic spacetime, which makes it very interesting for physics) For computer graphics, the 4dim "projective space" is very important. Here, rotations and translations fall together into one, more general group of operations. Instead of handling orientation and position separately and doing shifts and rotations in two separate steps, they can then be expressed very conveniently as one thing, in many cases even giving some performance gain. If you are interested, a good starting point is http://www.mrao.cam.ac.uk/~clifford/ but don't give up too quickly. At first it all seems rather strange, but after a short while, when you begin to understand, many things suddenly get a lot easier then they ever were in linear algebra.

Jul 04 2004

Norbert Nemec wrote:Stephen Waits wrote:I shall research the geometric algebra you mention. We generally just use quaternions for character animation, and possibly camera animation.

Geometric algebra just takes this to a more general level working for spaces of arbitrary dimension (also relativistic spacetime, which makes it very interesting for physics) For computer graphics, the 4dim "projective space" is very important. Here, rotations and translations fall together into one, more general group of operations. Instead of handling orientation and position separately and doing shifts and rotations in two separate steps, they can then be expressed very conveniently as one thing, in many cases even giving some performance gain.

I want to thank you for pointing me to this. We may be able to solve a problem we've run into with quaternions recently under this different algebra, using the reflection "tool".If you are interested, a good starting point is http://www.mrao.cam.ac.uk/~clifford/

I started with Jaap Suter's excellent tutorial located on his/her home page: http://jaapsuter.com/ There's quite a vocabulary with this new (to me) algebra; and I'll have to play with it quite a bit before I can claim any comfort - or even evaluate its usefulness in real-time 3D applications. There's a nice visualization program (GL based) which should help in my learning. Anyway, apologies to the group for going so OT. But, keep in mind, this is exactly the type of thing that we might want to do under D - cleanly and efficiently. Thanks again Norbert - you rock! --Steve

Jul 05 2004

Stephen Waits wrote:Norbert Nemec wrote:Geometric algebra...

I want to thank you for pointing me to this. We may be able to solve a problem we've run into with quaternions recently under this different algebra, using the reflection "tool".If you are interested, a good starting point is http://www.mrao.cam.ac.uk/~clifford/

I started with Jaap Suter's excellent tutorial located on his/her home page: http://jaapsuter.com/ There's quite a vocabulary with this new (to me) algebra; and I'll have to play with it quite a bit before I can claim any comfort - or even evaluate its usefulness in real-time 3D applications.

Please keep me informed on the outcome! I would really love to hear about real-world experience in that area.

Jul 06 2004

Take a look at my array suggestions at http://dlanguage.netunify.com/56 feel free to add to the wiki. Knud

Jul 02 2004

Knud Sørensen wrote:Take a look at my array suggestions at http://dlanguage.netunify.com/56 feel free to add to the wiki.

I think I've read these suggestions before. My feeling is that they are a bit too specialized to a certain kind of application, and that the simplicity of the syntax for the special cases you are listing would make it hard to grasp more general cases as well. Specifically: * The syntax you propose for "swizzling" and "write masking" collides badly with the multidimensional array syntax that has already started to become part of the language with the introduction of the new opIndex. * The Einstein summation in mathematics is just laziness. It makes sense to drop the summation sign in rotationally (or in spacetime: relativistic) invariant expressions. But outside of geometry, indices may mean many other thing than tensor indices. Turning everything into a sum just because there are two identical indices appearing would be really annoying there. I am thinking about a good syntax for index notation that would probably capture all the cases of use that you are mentioning. The result will probably need a bit more typing than your suggestions, but it should also be a lot more flexible.

Jul 03 2004

On Sat, 03 Jul 2004 10:06:59 +0200, Norbert Nemec wrote:Knud Sørensen wrote:Take a look at my array suggestions at http://dlanguage.netunify.com/56 feel free to add to the wiki.

I think I've read these suggestions before. My feeling is that they are a bit too specialized to a certain kind of application, and that the simplicity of the syntax for the special cases you are listing would make it hard to grasp more general cases as well. Specifically: * The syntax you propose for "swizzling" and "write masking" collides badly with the multidimensional array syntax that has already started to become part of the language with the introduction of the new opIndex. * The Einstein summation in mathematics is just laziness.

mathematically refactoring.It makes sense to drop the summation sign in rotationally (or in spacetime: relativistic) invariant expressions. But outside of geometry, indices may mean many other thing than tensor indices. Turning everything into a sum just because there are two identical indices appearing would be really annoying there.

You only make the summation when you make a multiplication so linear combinations of transformation is also possible. I have droped the co- and contra-variant notation as too domain specific so invariant transformation would have to be deal with explicit.I am thinking about a good syntax for index notation that would probably capture all the cases of use that you are mentioning. The result will probably need a bit more typing than your suggestions, but it should also be a lot more flexible.

Great, I look forward to it. Knud

Jul 03 2004

Knud Sørensen wrote:* The Einstein summation in mathematics is just laziness.

mathematically refactoring.

Of course - in analytic mathematics, it really hits the point. In a correct (i.e. invariant) expression, all the indices have to pair up exactly. Anyhow: linear algebra is just one use for indices in a computer language. There are many other uses as well where indices behave completely different.

Jul 03 2004

Why can't we make the current array operations serve this purpose?

Jul 02 2004

Walter wrote:Why can't we make the current array operations serve this purpose?

They are not flexible enough. Consider this example in pseudo-syntax: for(int i=0;i<10;i++) for(int j=0;i<10;i++) B[i,j] = sin(A[i,j]); How would you vectorize this? B[] = sin(A[]); would only work if sin is defined on arrays. So either you define a fixed set of functions in the language, or you define it in the library. But to define it in the library, you would still need to do the operation inside the function, and there you can only reside to a good-old-loop breaking all the performance gain of vectorization. One first idea might be to just use every function element-wise in array expressions, but then you would still run into problems on cases like: for(int i=0;i<10;i++) for(int j=0;i<10;i++) D[i,j] = A[i,j]+C[i]; What is needed is a "vector expression" that uses index notation. A few ideas, I'm shuffling around: index i,j; // declare i and j to be indices for a vector expression B[i,j] = sin(A[i,j]); D[i,j] = A[i,j] + C[i]; or B[] = [i=0..10,j=0..10](sin(A[i,j])); D[] = [i=0..10,j=0..10](A[i,j] + C[i]); In the second notation you could do really powerful things like: real[8,8] E = [i=0..8,j=0..8]((i+j)%2); // initialize a chess-board real[8,8] F = [i=0..8,j=0..8](sin(i*j)); // initialize a 2D-sin wave real[8,8] G = [i=0..8,j=0..8]((i+j)%2 ? sin(i*j) : 0); // initialize a masked sine-wave I'm not fully happy about the syntax yet. The ideas will probably just need some time to ripen. In any case, I believe the current array expressions can be preserved as well. Both kinds of expressions would merge together well into one concept. Just use the index-free notation for simple cases and use indices when you need the full flexibility.

Jul 03 2004

"Norbert Nemec" <Norbert.Nemec gmx.de> wrote in message news:cc5rm4$2juo$1 digitaldaemon.com...Walter wrote:Why can't we make the current array operations serve this purpose?

They are not flexible enough. Consider this example in pseudo-syntax: for(int i=0;i<10;i++) for(int j=0;i<10;i++) B[i,j] = sin(A[i,j]); How would you vectorize this? B[] = sin(A[]); would only work if sin is defined on arrays. So either you define a fixed set of functions in the language, or you define it in the library. But to define it in the library, you would still need to do the operation inside the function, and there you can only reside to a good-old-loop breaking

the performance gain of vectorization. One first idea might be to just use every function element-wise in array expressions, but then you would still run into problems on cases like: for(int i=0;i<10;i++) for(int j=0;i<10;i++) D[i,j] = A[i,j]+C[i]; What is needed is a "vector expression" that uses index notation. A few ideas, I'm shuffling around: index i,j; // declare i and j to be indices for a vector

B[i,j] = sin(A[i,j]); D[i,j] = A[i,j] + C[i]; or B[] = [i=0..10,j=0..10](sin(A[i,j])); D[] = [i=0..10,j=0..10](A[i,j] + C[i]); In the second notation you could do really powerful things like: real[8,8] E = [i=0..8,j=0..8]((i+j)%2); // initialize a

real[8,8] F = [i=0..8,j=0..8](sin(i*j)); // initialize a 2D-sin

real[8,8] G = [i=0..8,j=0..8]((i+j)%2 ? sin(i*j) : 0); // initialize a masked sine-wave I'm not fully happy about the syntax yet. The ideas will probably just

some time to ripen. In any case, I believe the current array expressions can be preserved as well. Both kinds of expressions would merge together well into one

Just use the index-free notation for simple cases and use indices when you need the full flexibility.

Sounds great! But the syntax really shouldn't look like that, the "[i=0..8,j=0..8]" part looks like an array litteral, even if they aren't declared like this. When i think about it again , it could work like that if array litterals had some other syntax. :)

Jul 03 2004

Ivan Senji wrote:Sounds great! But the syntax really shouldn't look like that, the "[i=0..8,j=0..8]" part looks like an array litteral, even if they aren't declared like this. When i think about it again , it could work like that if array litterals had some other syntax. :)

For one, there are no array literals at all. There are array initializers, but that's something different. Apart from that, it's just a first idea. The exact syntax of course has to be fit into the language much cleaner. In general, I'm not really sure, whether I like the idea of using brackets for both array literals and array indexing. Maybe, using braces for both array and struct literals might be a better idea. But that again might collide with block delimiters. It really is a pity that we have only three kinds of paired delimiters to choose from. :-(

Jul 03 2004

"Norbert Nemec" <Norbert.Nemec gmx.de> wrote in message news:cc5u5m$2nv8$1 digitaldaemon.com...Ivan Senji wrote:Sounds great! But the syntax really shouldn't look like that, the "[i=0..8,j=0..8]" part looks like an array litteral, even if they aren't declared like this. When i think about it again , it could work like that if array litterals had some other syntax. :)

For one, there are no array literals at all. There are array initializers, but that's something different.

Not now, but hopefully there will be in 2.0 :), they are mentioned in the future page.Apart from that, it's just a first idea. The exact syntax of course has to be fit into the language much cleaner. In general, I'm not really sure, whether I like the idea of using brackets for both array literals and array indexing. Maybe, using braces for both array and struct literals might be a better idea. But that again might collide with block delimiters. It really is a pity that we have only three kinds of paired delimiters to choose from. :-(

It really is a pity, but i don't think using [] for array literals and indexing is a bad thing as long as it is made unambiguous.

Jul 03 2004

"Norbert Nemec" <Norbert.Nemec gmx.de> wrote in message news:cc5rm4$2juo$1 digitaldaemon.com...Walter wrote: What is needed is a "vector expression" that uses index notation. A few ideas, I'm shuffling around: index i,j; // declare i and j to be indices for a vector

B[i,j] = sin(A[i,j]); D[i,j] = A[i,j] + C[i]; or B[] = [i=0..10,j=0..10](sin(A[i,j])); D[] = [i=0..10,j=0..10](A[i,j] + C[i]); In the second notation you could do really powerful things like: real[8,8] E = [i=0..8,j=0..8]((i+j)%2); // initialize a

real[8,8] F = [i=0..8,j=0..8](sin(i*j)); // initialize a 2D-sin

real[8,8] G = [i=0..8,j=0..8]((i+j)%2 ? sin(i*j) : 0); // initialize a masked sine-wave I'm not fully happy about the syntax yet. The ideas will probably just

some time to ripen.

Now that would be really cool - and useful :) I agree that the syntax has got to be thought about more though... but I cannot think of any alterate notation that would be as effective and compact as what you have proposed, so I reckon you're heading in the right direction, at least!

Jul 03 2004

Walter wrote:Why can't we make the current array operations serve this purpose?

We could do it - just not very optimally. For example, we CAN do this in C++; however, we have to jump through all kinds of hoops to get any efficiency at all - and even then there's a good chance you won't get all you need out of it. As it stands, D is no better than C++ in this regard; however, I believe it can be made considerably better without sacrifice. --Steve

Jul 05 2004

I think, you misunderstood the question by Walter: It is about the already specified "array operations" like: a[] = b[] + c[]; vs. the more general vector operations that I'm thinking about: a[] = [i=0..10,j=0..10](b[i,j] + c[i,j]); // syntax not decided yet!! The former already do offer more than C++ has: vectorizable code. The compiler can optimize very freely, since the order of execution is not specified. The latter do not offer much more, but they offer the same thing much more flexible. Especially, they offer the extension of vector expressions by using functions. Walter was just asking why the former kind of expressions does not suffice. I answered in a previous message. This situation is hardly comparable to C++, where you don't have any vectorizable expressions at all. Everything has to be done in the library via expression templates. Compared to that, D already is way ahead. It is just flexibility and extensibility that is lacking in the current solution. Stephen Waits wrote:Walter wrote:Why can't we make the current array operations serve this purpose?

We could do it - just not very optimally. For example, we CAN do this in C++; however, we have to jump through all kinds of hoops to get any efficiency at all - and even then there's a good chance you won't get all you need out of it. As it stands, D is no better than C++ in this regard; however, I believe it can be made considerably better without sacrifice.

Jul 05 2004

Norbert Nemec wrote:Hi there, just a few words about this nights thoughts. This is probably post-1.0, but it should already be considered now in order to avoid steps in the wrong direction. For high-performance numerics, vectorization of array operations is absolutely crucial. Currently, I know of several approaches to this problem in different languages: * HPF (High-Performance-Fortran) has powerful vector expression capabilities built in. AFAIK, these are always limited to one routine, i.e. you cannot spread vectorization across routines. This makes it mostly impossible to extend the existing set of vector operations. * Some very promising languages (e.g. Single-Assignment-C: http://www.sac-home.org) take the step to purely functional languages, or even purely declarative languages, where the compiler can optimize much more aggressively than in imperative languages.

I very much think this is the solution for the simple reason that it centers around giving the compiler more information, as opposed to pushing the burden onto programmers. Something that comes to mind is taking advantage of D's foreach and letting the optimizer be a bit more aggressive. The symbol 'foreach' already connotes the concept of doing a bunch of stuff in some unspecified order, so it's a good choice. For instance, given something like this, // not quite legal because of foreach float[4][] sourcePixels, destPixels; float[] alpha; foreach(float[4] src, inout float[4] dest, float a; sourcePixels, destPixels, alpha) { foreach(float s, inout float d; src, dest) { d = s * a + d * (1 - a); } } it's easy for both people and machines to see that, for the innermost loop, that no iteration depends on the result of any other iteration. So, hypothetically, the compiler could do things like break the expression up a bit, d *= (1 - a); d += s * a; then break those sub-expressions into two separate loops. From there, a decision to unroll or vectorize could be made. I think, though, that if it were this easy, it would have been done by now, so it might be best to ignore everything I just said. :) -- andy

Jul 03 2004

Well, I did think about using foreach for the concept of vectorization - anyhow: even though the specs are not very specific about this, foreach loops seem to have a defined order of execution. Just looking at the way opApply is implemented makes the whole thing more like a fancy notation for a good old loop. Trying to integrate this with the concept of vectorization seems more like mixing two different things into one. It is true that we should avoid new language concepts wherever existing ones can be extended, but it would be equally bad to try to force new meanings into existing concepts where they do not fit. Andy Friesen wrote:Norbert Nemec wrote:Hi there, just a few words about this nights thoughts. This is probably post-1.0, but it should already be considered now in order to avoid steps in the wrong direction. For high-performance numerics, vectorization of array operations is absolutely crucial. Currently, I know of several approaches to this problem in different languages: * HPF (High-Performance-Fortran) has powerful vector expression capabilities built in. AFAIK, these are always limited to one routine, i.e. you cannot spread vectorization across routines. This makes it mostly impossible to extend the existing set of vector operations. * Some very promising languages (e.g. Single-Assignment-C: http://www.sac-home.org) take the step to purely functional languages, or even purely declarative languages, where the compiler can optimize much more aggressively than in imperative languages.

I very much think this is the solution for the simple reason that it centers around giving the compiler more information, as opposed to pushing the burden onto programmers. Something that comes to mind is taking advantage of D's foreach and letting the optimizer be a bit more aggressive. The symbol 'foreach' already connotes the concept of doing a bunch of stuff in some unspecified order, so it's a good choice. For instance, given something like this, // not quite legal because of foreach float[4][] sourcePixels, destPixels; float[] alpha; foreach(float[4] src, inout float[4] dest, float a; sourcePixels, destPixels, alpha) { foreach(float s, inout float d; src, dest) { d = s * a + d * (1 - a); } } it's easy for both people and machines to see that, for the innermost loop, that no iteration depends on the result of any other iteration. So, hypothetically, the compiler could do things like break the expression up a bit, d *= (1 - a); d += s * a; then break those sub-expressions into two separate loops. From there, a decision to unroll or vectorize could be made. I think, though, that if it were this easy, it would have been done by now, so it might be best to ignore everything I just said. :) -- andy

Jul 03 2004