/* @TITLE "Matrix multiply"  */
/* 
 * matrix multiply example
 * sequential version
 *
 * David Kotz 1991
 */

#include <stdio.h>
#include "dfk.h"
#include "timing.h"

char *progname;

static void Initialize();
static void Multiply();
static void PrintMatrix();
char *alloc();

int N;					/* the size of each square matrix */
float **A, **B, **C;		/* the matrices */
boolean printmatrix;		/* should we print the matrix? */

/* @SUBTITLE "Main program: testing" */
main(argc, argv)
	int argc;
	char **argv;
{
    unsigned int start;		/* time in 100ths of seconds */
    double time;			/* elapsed time for multiply */

    /* Get arguments */
    progname = argv[0];

    if (argc < 2) {
	   fprintf(stderr, "usage: %s N [print]\n", progname);
	   exit(-1);
    }

    N = atoi(argv[1]);
    if (N <= 1) {
	   fprintf(stderr, "N=%d; should be > 1. (N=atoi(%s))\n", 
			 progname, N, argv[1]);
	   exit(-1);
    }

    printmatrix = (argc >= 3);

    /* INITIALIZE */

    (void) timer_init();

    A = (float **)alloc2d(N, N, sizeof(float));
    B = (float **)alloc2d(N, N, sizeof(float));
    C = (float **)alloc2d(N, N, sizeof(float));

    if (!(A && B && C)) {
	   fprintf(stderr, "%s: Not enough memory.\n", progname);
	   exit(1);
    }

    /* Initialize each matrix */
    Initialize();

    if (printmatrix) {
	   PrintMatrix("A", A, N, N);
	   PrintMatrix("B", B, N, N);
    }

    /* now compute the value of each result location */
    start = timer_get();

    Multiply();    /* N^3 flops */

    time = (timer_get() - start) * SECperTIC;
    
    if (printmatrix)
	 PrintMatrix("C", C, N, N);
    
    if (printmatrix)
	 printf("Time for multiply: %g sec  (%g flops) N=%d\n",
		   time, N*N*N/time, N);
    printf("1 %g %g %d\n", time, N*N*N/time, N);
}

/* @SUBTITLE "Initialize: initialize matrices" */
static void
Initialize()
{
    int row;				/* the row number */
    int col;				/* the column number */
    float *Arow;			/* a row of A */
    float *Brow;			/* a row of B */

    for (row = 0; row < N; row++) {
	   Arow = A[row];
	   Brow = B[row];
	   srandom((int)(row * (timer_get())));
	   for (col = 0; col < N; col++) {
		  Arow[col] = (random() % 10000 - 5000) / 50.;
		  Brow[col] = (random() % 10000 - 5000) / 50.;
	   }
    }
}

/* @SUBTITLE "Multiply: compute product" */
static void
Multiply()
{
    int row;				/* the row we are to compute */
    int col;				/* the column we are computing */
    int i;				/* the loop variable */
    double sum;			/* the running sum */
    float *Arow;			/* a row of A */

    /* this loop takes N*N*N flops */
    for (row = 0; row < N; row++) {
	   Arow = A[row];
	   
	   for (col = 0; col < N; col++) {
		  /* this loop takes N flops */
		  for (i = 0, sum = 0.0; i < N; i++) {
			 /* this is one flop */
			 sum += Arow[i] * B[i][col];
		  }
		  C[row][col] = sum;
	   }
    }
}

/* @SUBTITLE "PrintMatrix: print out a matrix" */
static void
PrintMatrix(name, A, m, n)
	char *name;
	float **A;			/* the matrix A[m][n] */
	int m,n;
{
    int i,j;

    printf("\nMatrix %s[%d][%d]:\n", name, m, n);
    for (i = 0; i < m; i++) {
	   printf("row %d\n", i);
	   for (j = 0; j < n; j++) {
		  printf("%g\t", A[i][j]);
	   }
	   printf("\n");
    }
}



/* @SUBTITLE "alloc" */
/* define alloc, used by alloc2d, to get memory */
char *
alloc(size)
	int size;
{
    return(malloc(size));
}