matrix mult.

[Daniel R. Levy]

I have some comments/questions regarding the ongoing discussion of matrix
multiplication and articles which have been posted thereon.  Those in
net.lang.f77 sit tight; f77 comes in near the end of the article.

In article <267 at h.cs.cmu.edu>, rfb at h.cs.cmu.edu (Rick Busdiecker) writes:
>
>char *valloc ();
>
>	[Code omitted]

Flames will come quickly if I am mistaken, but I can't seem to find any
valloc() on our Sys5R2 3B20.

In article <960 at turtlevax.UUCP>, ken at turtlevax.UUCP (Ken Turkowski) writes:
>In article <890 at ncoast.UUCP> simpsong at ncoast.UUCP (Gregory R. Simpson @ The North Coast) writes:
>>Below is a kludge version of a matrix multiplier in C. The primary
>>reason for posting this is in hope that someone out there knows the
>>proper method of doing it.
>Anyone who's diddled with two dimensional arrays, and has tried to think
>about how a compiler would do it has come up with something like the
>following:

[reduced from the original shar archive format]

>/* Matcat multiplies two (n x n) matrices together.  The source matrices
> * are given in A and B, and the result is returned in C.
> */
>
># ifndef ARRAYS
>
>matcat(C, A, B, n)
>register double *C;
>double *A, *B;
>register int n;
>{
>	register double *ap, *bp;
>	register int k, j, i;
>
>	for (i = n; i-- > 0; A += n) {		/* Each row in A */
>	    for (j = 0; j < n; j++) {		/* Each column in B */
>		ap = A;		/* Left of ith row of A */
>		bp = B + j;	/* Top of jth column of B */
>		*C = 0.0;
>		for (k = n; --k >= 0; bp += n)
>		    *C += *ap++ * (*bp);   /* *C += A[i'][k'] * B[k'][j]; */
		    ^^^^^

	What is this supposed to do with respect to the returned ARRAY?
	I don't ever see C getting bumped up to the next element.  Is
	this supposed to be *(C++) += ..... ?  I don't get it even then....
	(I don asbestos suit just in case)

>	    }
>	}
>}
>
># else
>
>matcat(C, A, B, n)
>register double *A, *B, *C;
>register int n;
>{
>	register int i, j, k;
>	double sum;
>
>	for (i = 0; i < n; i++) {
>	    for (j = 0; j < n; j++) {
>		sum = 0.0;
>		for (k = 0; k < n; k++)
>		    sum += A[n*i+k] * B[n*k+j];
>		C[n*i+j] = sum;
>	    }
>	}
>}
>
># endif

This last seems more reasonable (at least it works!).  In fact, it can be
easily hacked into a form which can be called on by f77 programs [and
which is probably more efficient that anything that could be written in
f77 itself due to the register variables and the lousy optimization :-(]:

/* call from f77 by:

	double precision a, b, c
	integer m
	data m /<whatever>/
C
C  Multiply m*m matrices a and b to arrive at matrix c.
C
	call matcat(c,a,b,m)
 */

matcat_(C, A, B, m)
register double *A, *B, *C;
int *m;
{
	register int n;
	register int i, j, k;
	double sum;
	n = *m;
	for (i = 0; i < n; i++) {
	    for (j = 0; j < n; j++) {
		sum = 0.0;
		for (k = 0; k < n; k++)
		    sum += A[n*i+k] * B[n*k+j];
		C[n*i+j] = sum;
	    }
	}
}

There are still some problems with this the way it is stated;
Fortran expects arrays to be in row-major order and C puts them
in column-major order (did I get it right, Guy?  You wrote an
article about this some time ago).  Also this does not deal with
non-square matrix multiplication.  It sounds like an easy enough
exercise to come up with such a f77-compatible matrix multiplication
routine, but does anyone out there have source handy already?