#include "header.h"

gsl_matrix * Create_Rand_Matrix(int m, int n){
	gsl_matrix *A;
	int i;
	//Allocate matrix
	A = gsl_matrix_alloc(m, n);
	//Fill with random numbers NOTE *A and alt form A->data[i]
	for(i=0; i<m*n; i++){
		(*A).data[i] = drand48();
		}
	return A;
	}
	
double gsl_dgemm_bench(int m, int n, int p, int runs){
	gsl_matrix *A, *B, *C;
	int 	 	i;
	clock_t  	t1, t2;
	double 	time = 0.0;
	
	A = Create_Rand_Matrix(m,p);
	B = Create_Rand_Matrix(p,n);
	C = Create_Rand_Matrix(m,n);
	
	t1 = clock();
	for(i=0; i<runs; i++){
     		gsl_blas_dgemm(CblasNoTrans, CblasNoTrans,
     	 		1.0, A, B, 0.0, C);     		
     		}
	t2 = clock();
	time= ((double)(t2-t1))/CLOCKS_PER_SEC;
	time = time/runs;
	
	return time;
	}

double myMult_bench(int m, int n, int p, int runs){
	gsl_matrix *A, *B, *C;
	int 	 	i;
	clock_t  	t1, t2;
	double 	time = 0.0;
	
	A = Create_Rand_Matrix(m,p);
	B = Create_Rand_Matrix(p,n);
	C = Create_Rand_Matrix(m,n);
	
	
	t1 = clock();
	for(i=0; i<runs; i++){
     		my_dgemm(A, B, C);     		
     		}
	t2 = clock();
	time= ((double)(t2-t1))/CLOCKS_PER_SEC;
	time = time/runs;
	
	return time;
	}
	
void 	my_dgemm(gsl_matrix * A, gsl_matrix * B, gsl_matrix * C){
	int i, j, k;
	int m=(*A).size1, n=(*B).size2, p=(*A).size2;
	double *ap=(*A).data,*bp=(*B).data,*cp=(*C).data;
	double sum;
	
	for(i=0;i<m;i++){
	 for(j=0;j<n;j++){
	  sum = 0.0;
	   for(k=0;k<p;k++){
	    sum += ap[i*p+k]*bp[k*n+j];  
                    // A[i][k]*B[k][j]  A is m by p and B is p by n 
	    cp[i*n+j] = sum;
                   }
	  }
	 }
	}