simpsong@ncoast.UUCP (Gregory R. Simpson @ The North Coast) (11/14/85)
*** Line eater poison ***
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. The problem I encountered was that in order
to pass a two dimensional array (hence a matrix) I needed to specify
the second bound of the array. However, I want this to be a generic
matrix multipler, for matrices of all sizes not just a fixed size. Since
the column increases most rapidly in a two dimensional array, declaring
a large array (like mat[][20]) puts all of a 3x3 into the first row.
In the much dreaded language of Fortran, one may simply pass a variable
to the dimension statement... this doesn't work in C. I believe this
is a job for pointers, but I was unsuccessful in getting them to work.
Please give me pointers on pointers (pun intended) for this application.
(or if you are ambitious send me a better version of code)
My fixes were fairly, crude, and the result is returned in a mxn vector
as opposed to an mxn matrix. Please send any pointers, references, code
help, flames to the address below. I have included my version below.
p.s. Once this issue is resolved, I will post an enhanced matrix
manipulation set, (mult., inversion, determinant, co-factors, etc.)
Thanks in advance,
Gregory R. Simpson
UUCP: decvax!cwruecmp!ncoast!simpsong
CSNET: simpsong@case.CSNET
ARPA: simpsong@purdue-ecn.ARPA
p.s. I am also looking for the manual page to roff.c posted from site
gloria about 2 months back...
======= bend - mutilate - cut - line-eat - here - ok - ================
/* a sample main program to call matmult */
main()
{
float a[3][3],b[3][3];
float c[9];
int i,j;
/* a typical array */
a[0][0]=1.0;
a[0][1]=2.0;
a[0][2]=3.0;
a[1][0]=4.0;
a[1][1]=5.0;
a[1][2]=6.0;
a[2][0]=7.0;
a[2][1]=8.0;
a[2][2]=9.0;
/* another typical array */
for (i=0;i<=2;i++) {
for (j=0;j<=2;j++) {
b[i][j] = a[i][j];
}
}
/* multiply two typical arrays together */
matmlt(a,b,c,3,3,3);
/* print the result, (it is in a 3 times 3 vector... boo hiss) */
for (i=0;i<=8;i++) {
printf("%f\n",c[i]);
}
}
/* the matrix multiplier */
matmlt(a,b,c,n,m,l)
/* this is a lot easier in fortran, cause then I can say dimension a(m,n) */
int n, m, l;
float a[], b[];
float c[];
{
float a2[10][10], b2[10][10], c2[10][10];
int i,j,k,inc;
/* load the arrays I passed into an array - kludgy ehhh? */
inc = 0;
for (i=0;i<=2;i++) {
for (j=0;j<=2;j++) {
a2[i][j] = a[inc];
b2[i][j] = b[inc];
inc++;
}
}
for (i = 0; i <= n-1 ; i++) {
for (k = 0; k <= l-1 ; k++) {
c2[i][k] = 0;
for (j = 0; j <= m-1 ; j++) {
c2[i][k] = (a2[i][j])*(b2[j][k]) + c2[i][k];
}
}
}
/* load the array result back into a vector to pass back */
inc = 0;
for (i=0;i<=2;i++) {
for (j=0;j<=2;j++) {
c[inc] = c2[i][j];
printf("%f %f %f \n",a2[i][j],b2[i][j],c2[i][j]);
inc++;
}
}
}rfb@h.cs.cmu.edu (Rick Busdiecker) (11/17/85)
Subject: Re: matrix mult. (source and question)
Message-Id: <501087498/rfb@H.CS.CMU.EDU>
In-Reply-To: Gregory R. Simpson @ The North Coast's netnews message of Thu, 14-Nov-85 00:18:16 EST
If Gregory Simpson's data representation were slightly different, things
would go much more nicely, maybe even more nicely than in fortran 8-}
Declaring a two-dimensional array is not equivalent to declaring an
array of one-dimensional arrays, i.e.
double matrix [m][n]
is not the same as
double (matrix [m]) [n];
however a reference to element (m,n) may still be coded as matrix[m][n] in
a routine regardless of which declaration is used. With the second
declaration element (m,n) could also be referenced as *(*(matrix + m) + n).
/*===========================================================================
* matrix_multiply
*
* Multiplies the matrix a (dimension m,n) by the matrix b (dimension n,p)
* into the matrix c (dimension m,p). Matrices must be declared explicitly
* to be arrays of arrays rather than 2d arrays.
*/
void matrix_multiply (a, b, c, m, n, p)
double **a, **b, **c;
int m, n, p;
{
int i, j, k;
for (i = 0; i < m; i++)
for (j = 0; j < p; j++)
{
c [i][j] = 0.0;
for (k = 1; k < n; k++)
c += a [i][k] * b [k][j];
}
}ken@turtlevax.UUCP (Ken Turkowski) (11/18/85)
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: ----------------------------------------------------------------- # This is a shell archive. Remove anything before this line, then # unpack it by saving it in a file and typing "sh file". (Files # unpacked will be owned by you and have default permissions.) # # This archive contains: # matcat.c echo x - matcat.c cat > "matcat.c" << '//E*O*F matcat.c//' /* 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]; */ } } } # 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 //E*O*F matcat.c// exit 0 -- Ken Turkowski @ CIMLINC (formerly CADLINC), Menlo Park, CA UUCP: {amd,decwrl,hplabs,seismo,spar}!turtlevax!ken ARPA: turtlevax!ken@DECWRL.DEC.COM
hester@ICSE.UCI.EDU (Jim Hester) (11/20/85)
Your method is metter in terms of space, and can usually be made better in terms of time. There are some minor improvements, and some pro's and cons to consider. The basic change your code made, relative to the original procedure, was to have the programmer explicitly deal with his own linear array representation. This is a bit more efficient space-wise (than the pointer-vector representation), makes creating a new array much easier, and can probably be made more efficient time-wise than a compiler since, although a compiler might do static array lookups a bit faster, your code could be modified to pre-calculate some terms outside of loops for moderate speedups, which only an extremely intelligent optimizing compiler (if any) could do. The major disadvantage of this scheme is that the programmer cannot use the facilities already provioded by the language for manipulating multi-dimensional arrays, which leads to slight conceptual confusion requiring a fair amount of documentation explaining to potential modifiers your method and why they must uniformly conform to it. The programmer must also implement and consistantly use procedures/macros for array access, changing all lines in all procedures which make use of the arrays from (the system's notation): i = A[ x ][ y ]; A[ x ][ y ] = i; to: i = lookup( A, n, x, y ); assign( A, n, x, y, i ); or: i = A[ loc(n, x, y) ]; A[ loc(n, x, y) ] = i; (I generally use the second method, since it only requires one function, is conceptually closer to what the compiler does, is visually closer to the language's notation, is simpler, and much more easily allows the use of macros.) The original code was intended to be as general as possible, so the program using explicit linear arrays and my examples above should be modified to support n by m arrays (trivial). If you want to limit the problem to square arrays and believe they could get larger than about 50 by 50, there is a reasonably non-complex, asymptotically faster algorithm. Strassen's multiplication algorithm (which I'll spare the readers: see any decent algorithms text) takes Order( n^2.81 ) time as opposed to Gaussian algorithms (all those presented recently), which take Order( n^3 ) time.
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?