levy@ttrdc.UUCP (Daniel R. Levy) (11/21/85)
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@h.cs.cmu.edu>, rfb@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@turtlevax.UUCP>, ken@turtlevax.UUCP (Ken Turkowski) writes: >In article <890@ncoast.UUCP> simpsong@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?