#ifdef TIME
#include <system.cf>	/* Get DOS clock function. */
#else
long clock() {return 0;}
#endif
#include <stdio.h>

extern void * malloc(unsigned int);
extern void free(void *);

#define Trace_malloc 1

unsigned Size_mem() {
   unsigned int Alloc_amount 
   	= sizeof(int) == 2 ? 60000 : 0xff_ffff /* 16 mb. */,
   	Sum = 0;
   void try() {
      void *p;
      do {
         p = malloc(Alloc_amount);
         Alloc_amount = Alloc_amount / 2;
      } while (p == 0 && Alloc_amount > 256);  
      if (Trace_malloc) printf("Malloc-ed %d: %x\n",Alloc_amount*2,p);   
      if (p != 0) {
      	 Sum += Alloc_amount*2;
         if (Alloc_amount > 256) try();
         }
      if (p != 0) free(p);
      }
   try();
   return Sum;
   }    

char YesNo(char *msg) {
   char response;
   printf("%s (y/n)?",msg);
   scanf(" %c",&response);
   return response == 'y' || response == 'Y';
   } 

void print_banner() {
   int i;
	printf("\t\311");for(i=0;i<58;i++) printf("\315");printf("\273\n");
	printf("\t\272                  Big Matrix Multiply                     \272\n");
	printf("\t\314");for(i=0;i<58;i++) printf("\315");printf("\271\n");
	printf("\t\272  This program multiplies a square matrix by itself.      \272\n");
	printf("\t\272  The matrix contains 32-bit integers.  Then the program  \272\n");
	printf("\t\272  does the same for a matrix containing 32-bit floating-  \272\n");
	printf("\t\272  point numbers. You can specify the size of the matrix.  \272\n");
	printf("\t\272  M[i,j] is initialized to i+j.                           \272\n");
	printf("\t\310");for(i=0;i<58;i++) printf("\315");printf("\274\n");
   }

void main (int argc, char** argv) {
   int Size;
   long int *Matrix_in, *Matrix_out;
   float *Fmatrix_in, *Fmatrix_out;
/*   long int (*M100_in)[100][100], (*M100_out)[100][100];  */
   void init() {
#     define DEFAULT 100
      /* Get the parameters. */
      /* If one or more arguments, assume default size. */
      if (argc > 1) Size = DEFAULT;
      else {
      L: printf("Enter the number of rows (and columns) in the matrix.\n"
      		"Enter 0 for the default of %d.\n"
                "Now enter your number:",DEFAULT,argv);
      	 scanf("%d",&Size);
      	 if (Size == 0) Size = DEFAULT;
      	 if (Size < 0) {
      	    printf("Sorry, please try something a little bigger.\n");
      	    goto L;
      	    }
         }
   L2:printf("\nAllocating two %d by %d matrices of 32-bit quantities...\n", 
   			Size,Size);
      Matrix_in = malloc(Size*Size*sizeof(long int));
      if (Trace_malloc) printf("Malloc-ed %d: %x\n",Size*Size,Matrix_in);
      Matrix_out = malloc(Size*Size*sizeof(long int));
      if (Trace_malloc) printf("Malloc-ed %d: %x\n",Size*Size,Matrix_out);
      if (Matrix_in == 0 || Matrix_out == 0) {
      	 extern double sqrt(double);
      	 if (Matrix_in != 0) free(Matrix_in);
      	 printf("Sorry, I don't seem to have that much memory...hang on...\n");
      	 printf("I seem to be able to allocate up to %d bytes for my array.\n",
      	 	Size = Size_mem());
      	 printf("Since I need two matrices, the maximum number of rows is %d.\n",
      	 	Size = sqrt((double)(Size/sizeof(long int)/2)));
      	 if (YesNo("Do you want me to use all of memory")) goto L2;
      	 else goto L;
      	 }
      }
#define Subscript(M,R,C) M[(R)*Size+(C)]
   int Row, Col;
   print_banner();   
   init();

   long start_time = clock();
   printf("\nInitializing matrix...\n");
   for (Row = 0; Row < Size; Row++)
      for (Col = 0; Col < Size; Col++)
         Subscript(Matrix_in,Row,Col) = Row+Col;
         
   printf("\nMultiplying matrix...\n");
   printf("\n(It will require %ld individual element multiplications.)\n"
   	  " and %ld multiplications just to compute array subscripts.)\n",
   		Size*Size*Size,Size*Size*Size*2+Size*Size);
   for (Row = 0; Row < Size; Row++)
      for (Col = 0; Col < Size; Col++) {
         /* Multiply M[Row,*] by M[*,Col] to get Result[Row,Col]. */
         long int Sum = 0; int R2;
         for (R2 = 0; R2 < Size; R2++) 
            Sum += Subscript(Matrix_in,Row,R2)*Subscript(Matrix_in,R2,Col);
         Subscript(Matrix_out,Row,Col) = Sum;
         }   

   if (start_time != 0) 
      printf("That took %d hundredths of a second.\n\n",
   		clock()-start_time);
   start_time = clock(); 
   		
   printf("\nNow for floating point.  Initializing matrix...\n");
   Fmatrix_in = (float *) Matrix_in;
   Fmatrix_out = (float *) Matrix_out;
   for (Row = 0; Row < Size; Row++)
      for (Col = 0; Col < Size; Col++)
         Subscript(Fmatrix_in,Row,Col) = Row+Col;
         
   printf("\nMultiplying matrix...\n");
   for (Row = 0; Row < Size; Row++)
      for (Col = 0; Col < Size; Col++) {
         /* Multiply M[Row,*] by M[*,Col] to get Result[Row,Col]. */
         float Sum = 0; int R2;
         for (R2 = 0; R2 < Size; R2++) 
            Sum += Subscript(Fmatrix_in,Row,R2)*Subscript(Fmatrix_in,R2,Col);
         Subscript(Fmatrix_out,Row,Col) = Sum;
         }   
   if (start_time != 0) 
   printf("That took %d hundredths of a second.\n\n",
   		clock()-start_time);
   
   printf("\nEnd of matrix multiply.\n");
   }
