Relay-Version: version B 2.10 5/3/83; site utzoo.UUCP
Posting-Version: version B 2.10.1 6/24/83; site sdcsla.UUCP
Path: utzoo!watmath!clyde!burl!mgnetp!ihnp4!zehntel!dual!amd70!decwrl!decvax!ittvax!dcdwest!sdcsvax!sdcsla!west
From: west@sdcsla.UUCP
Newsgroups: net.physics
Subject: Re: VMS vs. Unix (programs)
Message-ID: <590@sdcsla.UUCP>
Date: Tue, 3-Jul-84 00:03:33 EDT
Article-I.D.: sdcsla.590
Posted: Tue Jul  3 00:03:33 1984
Date-Received: Sat, 30-Jun-84 03:12:22 EDT
References: <360@noao.UUCP>
Organization: UC San Diego: Institute for Cognitive Science
Lines: 107

To make things a little fairer, I rewrote the C matrix multiplication
program to take advantage of C's pointers.   Of course, this is something a
good ``optimizer'' would do for you, but it isn't that hard to do by hand.

Peter Allan's original program:

		/*  Multiply two matrices  */

		#define N 100

		main()
		{
		   float a[N][N],b[N][N],c[N][N],sum;
		   int i,j,k;

		   /*  Ininialise a and b  */

		   for (i=0 ; i<N ; i++)
		      for (j=0 ; j<N ; j++)
		      {
			 a[i][j] = i+j;
			 b[i][j] = i-j;
		      }

		/*  Multiply a by b to give c  */

		   for (i=0 ; i<N ; i++)
		      for (j=0 ; j<N ; j++)
		      {
			 sum = 0.0;
			 for (k=0 ; k<N ; k++)
			    sum += a[i][k]*b[k][j];
			 c[i][j] = sum;
		      }
		}


My quick and easy improvements:

	/*  Multiply two matrices  */

	#define N 100

	main()
	{
	   float a[N][N],b[N][N],c[N][N];
	   register float sum;
	   register float *pa, *pb, *pc;
	   register int i,j,k;

	   /*  Initialize a and b  */

	   for (i=0 ; i<N ; i++)
	      {
	      pa = a[i];		/* actually could be done just once */
	      pb = b[i];
	      for (j=0 ; j<N ; j++, pa++, pb++)
		  {
		  *pa = i+j;
		  *pb = i-j;
		  }
	      }

	/*  Multiply a by b to give c  */

	   for (i=0 ; i<N ; i++)
	      {
	      pa = a[i];
	      pb = b[i];
	      for (j=0 ; j<N ; j++, pa++, pb++)
		 {
		 sum = 0.0;
		 for (k=0 ; k<N ; k++)
		    sum += (*pa)*(*pb);
		 c[i][j] = sum;
		 }
	      }
	}

---------------

So the differences are simply judicious use of "register" declarations
and using pointers and pointer incrementation to step through an array.
The initializations of "pa" and "pb" could even be moved out of the
outer loops, but for very little gain.

Timings taken on a fairly uncrowded VAX 750 (no floating point extras),
for ``optimized'' compilation:

time mm1	-- original C program, "optimized".
67.8u 11.6s 2:17 57% 1+202k 0+0io 3pf+0w
67.3u 6.2s 2:31 48% 1+201k 2+1io 1pf+0w

time mm1a	-- improved C program, "optimized".
21.9u 2.5s 0:40 60% 1+198k 0+0io 2pf+0w
21.9u 2.2s 0:40 60% 1+200k 0+0io 2pf+0w

time mm2	-- original Fortran program, "optimized".
28.1u 3.5s 0:51 61% 11+243k 1+1io 11pf+0w
28.6u 3.0s 1:09 45% 11+244k 0+0io 5pf+0w

Note that now C at least beats f77 on its (C's) home turf.

		-- Larry West, UC San Diego, Institute for Cognitive Science
		-- decvax!ittvax!dcdwest!sdcsvax!sdcsla!west
		-- ucbvax!sdcsvax!sdcsla!west
		-- west@NPRDC