/****Copyright (c) 2001  Alexandros V. Gerbessiotis. 
 **** All rights reserved.                               
 ****Permission to use, copy, modify, and distribute this software,
 ****and to incorporate it, in whole or in part, into 
 ****other software,is hereby granted without fee, provided that
 ****  (1) the above copyright notice and this permission 
 ****      notice appear in all copies of the source code, and 
 ****      the above copyright notice appear in clearly visible 
 ****      form on all supporting documentation and distribution 
 ****      media, and
 ****     (2) any redistribution of the software, in original
 ****         or modified form, be without fee and subject to
 ****         these same conditions.
 ****   No guarantee, whatsoever, is offered that the code works as
 ****   advertised or is absent of any, even damaging,
 ****   side-effects. Use this code at your own personal risk.
 ****/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <limits.h>
#define	TEN 8
#define MATDEBUG
#define MATMULONLY

extern clock_t clock();

void initialize(double **m4,int n,int debug,int runs);
void print_mat(double *mat,int n,int debug);
void mult(double *,double *,double *,int );
void mult2(double *,double *,double *,int );
void mult3(double *,double*,double*,int);
void fast();
void fast2();
void transpose();

int 	main(int argc, char *argv[])
{
int 	i,j,k,n,runs,debug;
double  *mat,*mat1,*mat2,*mat3,*temp;
clock_t t1,t3,t4;
double 	a,b,c;

 n=16;debug=0; runs=1;
 if (1 == argc ) {
  printf("Usage: %s  dimension runs debug \n",argv[0]);
  return(1);
 } 
 else  
  if (4 == argc){
    n=atoi(argv[1]); 
    runs =atoi(argv[2]);
    debug=atoi(argv[3]);
  } 

 for(i=0;i<argc;i++) printf("Parameters: %s ",argv[i]); printf("\n");

 initialize(&mat,n,debug,runs);
 print_mat(mat,n,debug);
 free((void *)mat);
}

void check_if_null(void *x,char *s,char *t)
{
  if (NULL == x) {
     printf("%s:: Null pointer error (%s)\n",s,t);fflush(stdout);
     abort();
  }
}

void initialize(double **mat_r,int n, int debug,int runs)
{
double *mat1,*mat2,*mat;
register int i,j,k;
double t1,t3,t4,a,b,c;

 mat1 = (double *) malloc(n*n*sizeof(double));
#ifdef MATDEBUG
 check_if_null(mat1,"initialize","mat1");
#endif
 mat2 = (double *) malloc(n*n*sizeof(double));
#ifdef MATDEBUG
 check_if_null(mat2,"initialize","mat2");
#endif

 for(j=0;j<n;j++)
   for(i=0;i<n;i++) {
    mat1[j*n+i] = (double) (j+1);
    mat2[j*n+i] = (double) (2*j+i+1);
   }
 if (1==debug){
  printf("Original Matrix A \n");
  for(i=0;i<n;i++)
    for(j=0;j<n;j++)
      printf("%6.2f%c",mat1[i+j*n],((j<(n-1))?' ':'\n'));
  printf("Original Matrix B \n");
  for(i=0;i<n;i++)
    for(j=0;j<n;j++)
      printf("%6.2f%c",mat2[i+j*n],((j<(n-1))?' ':'\n'));
 }


 mat = (double *) malloc(n*n*sizeof(double));
#ifdef MATDEBUG
 check_if_null(mat,"initialize","mat");
#endif
 for(j=0;j<n;j++)
  for(i=0;i<n;i++)
    mat[j*n+i]=0.0;

 t1=0.0;
 t3=clock();
   transpose(mat1,mat1,n);
   for(i=0;i<runs;i++) {
     mult(mat1,mat2,mat,n);
   }
 t4=clock(); 
 t1 += (t4-t3)/((double)CLOCKS_PER_SEC);

   a = (double)n*n*n ; 
   b = (double)n*n*n + (double) n*n*(n-1) ; 
   c = b;
 if (t1 != 0.0) {
  printf("Matrix multiplicatio1: %8.5f/%8.5f (Mflop: %8.3f) \n", t1/runs,t1, (runs*c)/(t1*1000000.0)); 
 } 
 else {
  printf("Matrix multiplicatio1:No time reported. Run too short\n");
 }

 if (1==debug){
   printf("Resultant Matrix of multiplication1\n");
   for(i=0;i<n;i++)
     for(j=0;j<n;j++)
       printf("%6.2f%c",mat[i+j*n],((j<(n-1))?' ':'\n'));
 }

 for(j=0;j<n;j++)
   for(i=0;i<n;i++) {
    mat1[j*n+i] = (double) (j+1);
    mat2[j*n+i] = (double) (2*j+i+1);
   }
 if (1==debug){
  printf("Original Matrix A \n");
  for(i=0;i<n;i++)
    for(j=0;j<n;j++)
      printf("%6.2f%c",mat1[i+j*n],((j<(n-1))?' ':'\n'));
  printf("Original Matrix B \n");
  for(i=0;i<n;i++)
    for(j=0;j<n;j++)
      printf("%6.2f%c",mat2[i+j*n],((j<(n-1))?' ':'\n'));
 }

 for(j=0;j<n;j++)
  for(i=0;i<n;i++)
    mat[j*n+i]=0.0;

 t1=0.0;
 t3=clock();
   transpose(mat1,mat1,n);
   for(i=0;i<runs;i++)
     mult2(mat1,mat2,mat,n);
 t4=clock(); 
 t1 += (t4-t3)/((double)CLOCKS_PER_SEC);

   a = (double)n*n*n ; 
   b = (double)n*n*n + (double) n*n*(n-1) ; 
   c = b;
 if (t1 != 0.0) {
  printf("Matrix multiplicatio2: %8.5f/%8.5f (Mflop: %8.3f) \n", t1/runs,t1, (runs*c)/(t1*1000000.0)); 
 } 
 else {
  printf("Matrix multiplicatio2:No time reported. Run too short\n");
 }

 if (1==debug){
   printf("Resultant Matrix of multiplication2\n");
   for(i=0;i<n;i++)
     for(j=0;j<n;j++)
       printf("%6.2f%c",mat[i+j*n],((j<(n-1))?' ':'\n'));
 }

 for(j=0;j<n;j++)
   for(i=0;i<n;i++) {
    mat1[j*n+i] = (double) (j+1);
    mat2[j*n+i] = (double) (2*j+i+1);
   }
 if (1==debug){
  printf("Original Matrix A \n");
  for(i=0;i<n;i++)
    for(j=0;j<n;j++)
      printf("%6.2f%c",mat1[i+j*n],((j<(n-1))?' ':'\n'));
  printf("Original Matrix B \n");
  for(i=0;i<n;i++)
    for(j=0;j<n;j++)
      printf("%6.2f%c",mat2[i+j*n],((j<(n-1))?' ':'\n'));
 }
 for(j=0;j<n;j++)
  for(i=0;i<n;i++)
    mat[j*n+i]=0.0;


 t1=0.0;
 t3=clock();
   transpose(mat1,mat1,n);
   for(i=0;i<runs;i++)
     mult3(mat1,mat2,mat,n);
 t4=clock(); 
 t1 += (t4-t3)/((double)CLOCKS_PER_SEC);

   a = (double)n*n*n ; 
   b = (double)n*n*n + (double) n*n*(n-1) ; 
   c = b;
 if (t1 != 0.0) {
  printf("Matrix multiplicatio3: %8.5f/%8.5f (Mflop: %8.3f) \n", t1/runs,t1, (runs*c)/(t1*1000000.0)); 
 } 
 else {
  printf("Matrix multiplicatio3:No time reported. Run too short\n");
 }

 if (1==debug){
   printf("Resultant Matrix of multiplication3\n");
   for(i=0;i<n;i++)
     for(j=0;j<n;j++)
       printf("%6.2f%c",mat[i+j*n],((j<(n-1))?' ':'\n'));
 }

#ifndef MATMULONLY
  t1=0.0;
  t3=clock();
   for(i=0;i<runs;i++)
    fast(mat,n);
  t4=clock(); 
  t1 += (t4-t3)/((double)CLOCKS_PER_SEC);

  b = (double)n*n*n *0.5 ; c=2.0*(n*n*n-n);
  if (t1 != 0.0){
    printf("Loop1 execution      : %8.5f/%8.5f (Mflop: %8.3f)[Ignore:%8.2f]\n", t1/runs,t1,
              (c)/(t1*1000000.0), (log(log(mat[0])))); 
  }
 else {
  printf("Loop1 execution      :No time reported. Run too short\n");
 }


  t1=0.0;
  t3=clock();
   for(i=0;i<runs;i++)
    mult2(mat1,mat2,mat,n);
  t4=clock(); 
  t1 += (t4-t3)/((double)CLOCKS_PER_SEC);
  free((void *) mat1); 
  free((void *) mat2); 

   a = (double)n*n*n ; 
   b = (double)n*n*n + (double) n*n*(n-1) ; 
   c = b;
  if (t1 != 0.0){
    printf("Loop2 execution      : %8.5f/%8.5f (Mflop: %8.3f)\n",t1/runs,t1,(c)/(t1*1000000.0)); 
  }
 else {
  printf("Loop2 execution      :No time reported. Run too short\n");
 }

#endif

  *mat_r = mat;
}

void	print_mat(double *mat,int n, int debug)
{
 register int i,j;

   if (1==debug){
     printf("Resultant Matrix \n");
     for(i=0;i<n;i++)
       for(j=0;j<n;j++)
         printf("%6.2f%c",mat[i+j*n],((j<(n-1))?' ':'\n'));
   }
}

