/*
Written by Holguer A Becerra
UBC Student
Master of Biomedical Engineering 
Course: ECE 328 
Prof. Guy Lemiuex.
Year 2015-Winter Term 1
*/


#include <stdio.h>
#include <sys/time.h>
#include <time.h>
#include <stdlib.h>
#include <string.h>
 

#define INT 0 
#define FLOAT 1 
#define DOUBLE 1 
#define MIN(a,b) (((a)<(b))?(a):(b))
#define TYPE_VAR FLOAT

#if TYPE_VAR == INT
	int * matrix1;
	int * matrix2;
	int * matrix3;
	int sum;
#elif TYPE_VAR == FLOAT
	float * matrix1;
	float * matrix2;
	float * matrix3;
	float sum;
#elif TYPE_VAR == DOUBLE
	double * matrix1;
	double * matrix2;
	double * matrix3;
	double sum;
#endif 


typedef unsigned long long int64;
unsigned measurement=0,ut,mt;
void init_timer();
static __inline__ int64 getticks(void);


int main(int argc, char *argv[])
{
	int64 tick,tick1;
	int N,rows,columns;
	int c,d,k;
	int print=0;
	int B=64;
	//initializing timer variables and timer
	init_timer();
	
	if(argc<=1){
		printf("too few arguments");
		return 0;
	}
	
	// the first argument has the number of rows, in this case the number of columns are the same as the number of rows.
	sscanf(argv[1],"%d",&N);
	sscanf(argv[2],"%d",&print);
	sscanf(argv[3],"%d",&B);

	/* initialize random seed: */
	srand (100);
	
	
	/*Get Timer*/
	tick = getticks();
	
	
  // Initialize Matrices	
	#if TYPE_VAR == INT
		matrix1=(int *)malloc((N+1)*N*sizeof(int));
	#elif TYPE_VAR == FLOAT
		matrix1=(float *)malloc((N+1)*N*sizeof(float));
	#elif TYPE_VAR == DOUBLE
		matrix1=(double *)malloc((N+1)*N*sizeof(double));
	#endif 

	printf("Initializing matrices\n");
	for (rows= 0; rows < (N+1)*N; rows++){
		#if TYPE_VAR == INT
			*(matrix1+rows)=rand() % 1000 + 1;			
			//printf("1[%d,%d]=%d\n",rows/N,rows,*(matrix1+rows));
			//printf("2[%d,%d]=%d\n",rows/N,rows,*(matrix2+rows));
		#elif TYPE_VAR == FLOAT
			*(matrix1+rows)=rand() % 10 + 1.5;
			//printf("1[%d,%d]=%f\n",rows/N,rows,*(matrix1+rows));
			//printf("2[%d,%d]=%f\n",rows/N,rows,*(matrix2+rows));
		#elif TYPE_VAR == DOUBLE
			*(matrix1+rows)=rand() % 1000 + 1.5;
			//printf("1[%d,%d]=%f\n",rows/N,rows,*(matrix1+rows));
			//printf("2[%d,%d]=%f\n",rows/N,rows,*(matrix2+rows));
		#endif 
		
	}

  
  // Performance  Finish and Measure
  /*Get Timer*/
  tick1 = getticks();
  measurement = (unsigned)((tick1-tick)/ut);
  
  printf("\n Initialization time: %d us\n",measurement);
  printf("Initialization Done, moving over Forward Elminination\n");
 
 
 if(print){
	 for (rows= 0; rows < N; rows++){
		for (d= 0; d < N+1; d++){
			printf("%0.2f	",*(matrix1+rows*(N+1)+d));
		}
		printf("\n");
	  }
 }
 
 
 
 /*Get Timer*/
	tick = getticks();
 
 
 
   // Gauss Forward Elminination
    float pivot=0.0;
	float pivot_n=0.0;
	
	int section=0;
	int section2=0;
	int count=0;
	for(c=0;c<N-1; c++){// pivot pointer		
	   pivot_n=(*(matrix1 + (c)*(N+1) +c));
			for(k=c+1; k<N; k++){
				pivot=(*(matrix1 + (k)*(N+1) +c))/pivot_n;
				count=0;
				section=0;
				while(section<N+1){
					for(d=section; d<MIN(section+B,N+1); d++){
						*(matrix1 + (k)*(N+1) +d)=*(matrix1 + (k)*(N+1) +d)-pivot*(*(matrix1 + (c)*(N+1) +d));
						count++;
					}
					section+=count;
				}
			}
	}
	
  // Performance  Finish and Measure
  /*Get Timer*/
  tick1 = getticks();
  measurement = (unsigned)((tick1-tick)/ut);
  
  printf("\n Forward Elminination time: %d us\n\n",measurement);
  
  if(print){
	 for (rows= 0; rows < N; rows++){
		for (d= 0; d < N+1; d++){
			printf("%0.2f	",*(matrix1+rows*(N+1)+d));
		}
		printf("\n");
	  }
  }
  
  /*Get Timer*/
  tick = getticks();
  // back  sustitution
  for(c=N-1; c>=0; c--){
	for(k=c+1; k<N; k++){
		(*(matrix1 + (c)*(N+1) +N))=(*(matrix1 + (c)*(N+1) +N))-((*(matrix1 + (c)*(N+1) +k))*(*(matrix1 + (k)*(N+1) +N)));
		
	}
	(*(matrix1 + (c)*(N+1) +N))=(*(matrix1 + (c)*(N+1) +N))/(*(matrix1 + (c)*(N+1) +c));
  }
  
  /*Get Timer*/
  tick1 = getticks();
  measurement = (unsigned)((tick1-tick)/ut);
  
  printf("\n Back Sustitution time: %d us\n",measurement);
  printf("\n\n");
  
  if(print){
	 for (rows= 0; rows < N; rows++){
		for (d= 0; d < N+1; d++){
			printf("%0.2f	",*(matrix1+rows*(N+1)+d));
		}
		printf("\n");
	  }
  }
	
  // free memory
  
  free(matrix1);
  
  


 
  return 0;
}

static __inline__ int64 getticks(void)
{
     unsigned a, d;
     asm volatile("rdtsc" : "=a" (a), "=d" (d));
     return (((int64)a) | (((int64)d) << 32));
}


void init_timer(){
     int64 tick,tick1;
     

     // ut is the divisor to give microseconds
     // mt gives milliseconds

     FILE *pf;
     int i,r,l,n=0;
     char s[100];

     // time how long it takes to get the divisors, as a test 
     tick = getticks();

     
     pf = fopen("/proc/cpuinfo","r");
     do {
      r=fscanf(pf,"%s",&s[0]);
      if (r<0) {
       n=5; break;
      } else if (n==0) {
       if (strcmp("MHz",s)==0) n=1;
      } else if (n==1) {
       if (strcmp(":",s)==0) n=2;
      } else if (n==2) {
       n=3;
      };
     } while (n<3);
     fclose(pf);

     l=strlen(s);
     s[l-4]=s[l-3];
     s[l-3]=s[l-2];
     s[l-2]=s[l-1];
     s[l-1]=(char)0;

     mt=atoi(s);
     s[l-4]=(char)0;
     ut=atoi(s);

     //printf("%s Mhz - ut = %u, mt = %u // hardcode these for your a CPU-specific binary ;-)\n",s,ut,mt);

}