Array extensions (again)

[queueing-rmail]

I think that I promised to post an example of how to do a "matrix multiply"
(as that term is defined in the numerical computing world) using the
array extensions that I have suggested for consideration by NCEG.

It's a tricky problem.

At first I though that this problem might provide a justification for
proposing one more new operator (along with the anti-decay operator I
already proposed) but after thinking about it a bit more, I convinced
myself that it would be enormously difficult to design an operator which
would do anything particularly useful in the context of the "matrix
multiply" problem, and that even if I could dream up a new primitive
(and generally useful) operator to help implement a matrix multiply,
the semantics of it would be so complicated as to be downright confusing.

So I think I'll just stick with the generic tools I've already got.

I will however add one small refinement to what I already suggested as
regards to extending the semantics of existing C operators (in a natural
and intutive way) to "objects" of arbitrary dimensionality.  Previously,
I suggested extensing the semantics of existing operators (in a simple
way) so that (for example) if you used `+' between two operands which
were "shape compatible" arrays of some arbitrary dimensionality, then
the result would be itself be an array whose "shape" would be compatible
with that of the two input operands.  I also suggested (although perhaps
not very clearly) that for two array type values to be "shape compatible"
they must have the same dimensionality.

I now suggest that for all operators like a= (where a is one of +, -, *,
/, %, ^, |, %, ~, <<, >>) the shape compatability rules should be relaxed
somewhat and that it should be permissible for the left operand to be
shape compatible with the *element type* of the right operand.  In such
cases the semantics of the operation would be that the "basic" operation
is carried out repeatedly, once for each paring of the left operand with
an element of the right operand.  Thus, given declarations like:

		double array[100];
		double sum;

we would be permitted to say:

		sum = 0.0;
		sum += array;

thereby computing the sum of all elements of `array' into `sum'.

Now, with this small refinement in mind, I'll provide the code to do a
matrix multiply using the anti-decay operator, and the other array
extensions for C which I have previously suggested.  For the sake of
simplicity (of the example) I'll just write the code to multiply two
fixed size 2-dimensional arrays.

	#define SIZE 10
	
	float input1[SIZE][SIZE];
	float input2[SIZE][SIZE];
	float output[SIZE][SIZE];
	
	void matrix_multiply ()
	{
		register unsigned int i, j;
		float temp_vector[SIZE];
	
		for (i = 0; i < SIZE; i++) {

			/* Collect the i'th column of `input1' and put it
			   into `temp_vector'.  */

			for (j = 0; j < SIZE; j++)
				temp_vector[j] = input1[i][j];

			for (j = 0; j < SIZE; j++) {

				output[i][j] = 0.0;

				/* Multiply the i'th column of `input1' by the
				   j'th row of `input2', sum the resulting
				   vector, then save the result.  */

				output[i][j] += temp_vector[] * input2[j][];
			}
		}
	}

Now please don't get on my case if this doesn't quite implement a proper
matrix multiply (as the numerical analysts understand it)  I'm relatively
ignorant of such things, so I may not have quite gotten the semantics
exactly right as far as the purists are concerned, but this example should
be sufficient to illustrate the idea anyway.

I feel quite sure that some people will look at this example, understand
it completely, and that claim (perhaps rightly) that there is no way this
code can be maximally optimal for their own pet vector machines.  I'll grant
in advance that such assertions are probably true for this particular vector
machine or that particular vector machine, but if the NCEG is looking for a
*generalized* set of extensions for C which may prove useful on a wide
variety of different vector machines, I doubt that it is possible to do much
better than this... at least not without breaking the existing rules of
ANSI C in some far more violent way.

Comments are welcomed.  Separately, if anyone has any trouble interpreting
this example, send me E-mail and I'll be glad to explain it in more detail.


// Ron ("Loose Cannon") Guilmette    uucp: ...uunet!lupine!segfault!rfg
//
// 	"On the one hand I knew that programs could have a compelling
// 	 and deep logical beauty, on the other hand I was forced to
// 	 admit that most programs are presented in a way fit for
// 	 mechanical execution, but even if of any beauty at all,
// 	 totally unfit for human appreciation."
// 						-- Edsger W. Dijkstra