#include <stdlib.h>
#include <stdio.h> 
#include <time.h>

#define get_ticks(var) {					   \
      unsigned int __a, __d;					   \
      asm volatile("rdtsc" : "=a" (__a), "=d" (__d));		   \
      var = ((unsigned long) __a) | (((unsigned long) __d) << 32); \
   } while(0)


int dgemm_(char *, char *, int *, int *, int *, double *,
	   double *, int *, double *, int *, double *, double *, int *);



int main( int argc, char *argv[] ){

   /* 
      // part 1 

   double a = 3.1415;
   int i, j;
   unsigned long ticksB4, ticksAFT, ticks; 

   get_ticks(ticksB4);   
      
   for( i=0; i<10000; i++ )
      for( j=0; j<10000; j++ )
	 a = a * 1.000000001;

   get_ticks(ticksAFT); 
   ticks = ticksAFT - ticksB4; 
   printf(" loop time= %.5ld\n", ticks );

   printf( "result = %.16g\n", a );
   */


   /*
   // part 2

   //   double a = 3.1415, b = 2.72818, c;
   //   double a = 3.1415, b, c;
   double a = 0, b, c;

   b = atof(argv[1]);

   c = a * b;

   printf( "result = %.16g\n", c );
   */



   
   // part 3: GEMM

   int m,n,k;
   int i,j,l;
   double *A, *B, *C, *D;
   unsigned long ticksB4, ticksAFT, ticks2; //ticks1; 

   m = atof(argv[1]);
   n = atof(argv[2]);
   k = atof(argv[3]);

   A = (double *) malloc( m*k * sizeof(double) );
   B = (double *) malloc( k*n * sizeof(double) );
   C = (double *) malloc( m*n * sizeof(double) );
   D = (double *) malloc( m*n * sizeof(double) );

   srand48( (unsigned)time((time_t *)NULL) ); 
   
   for( i=0; i<m; i++ )
      for( j=0; j<k; j++ ) 
	 A[i+j*m] = drand48();  

   for( i=0; i<k; i++ )
      for( j=0; j<n; j++ ) 
	 B[i+j*k] = drand48();  

   for( i=0; i<m; i++ )
      for( j=0; j<n; j++ ){
	 C[i+j*m] = 0;
	 D[i+j*m] = 0;
      }
      
   /*
   get_ticks(ticksB4);   
      
   for( i=0; i<m; i++ )
      for( j=0; j<n; j++ )
	 for( l=0; l<k; l++ )
	    C[i+j*m] += A[i+l*m] * B[l+j*k]; 


   get_ticks(ticksAFT); 
   ticks1 = ticksAFT - ticksB4; 
   printf(" triple loop: %.5ld\n", ticks1 );
   */

   int ITER = 100;
   
   for( i=0; i<ITER; i++ ){
      
      get_ticks(ticksB4);   
      
      double one = 1.0, zero = 0.0;
      
      dgemm_( "N", "N", &m, &n, &k, &one, A, &m, B, &k, &zero, D, &m ); 
      
      get_ticks(ticksAFT); 
      ticks2 = ticksAFT - ticksB4;
      printf(" iter(%d): %.5g\n", i+1, (double) 2*m*n*k / (double) ticks2 );
      //      printf(" BLAS: %.5ld\t\t\t ratio:%.3f\n", ticks2, (double) ticks1/(double) ticks2 );
   }
      
   
   /*
   for( i=0; i<m; i++ ){
      for( j=0; j<n; j++ )
	 printf( " %.16g,",  (C[i+j*m]-D[i+j*m])/D[i+j*m] );
      printf(" \n");
   }
   */
   
   
   return(0);
}