// LAPACK: -L/Users/pauldj/works/libs/lapack-3.5.0/  -llapack  -lgfortran
// BLAS:   -L/opt/local/lib -lopenblas

#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 dpotf2_(char *, int *, double *, int *, int *);
int dpotrf_(char *, int *, double *, int *, int *);

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

void mat_print( double *A, int n, int m, int ldA, char *text);



int main( int argc, char *argv[] ){
   int i, j, m, info;
   double *A, *A0;
   unsigned long ticksB4, ticksAFT, ticks; 

   m = atof(argv[1]);
   
   A = (double *) malloc( m*m * sizeof(double) );
   A0= (double *) malloc( m*m * sizeof(double) );

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

   //mat_print( A, m, m, m, "A" );
   //mat_print( A0, m, m, m, "A0" );


   int ITER = 10;
   
   for( i=0; i<ITER; i++ ){
      get_ticks(ticksB4);   
   
      //dpotf2_( "L", &m, A, &m, &info );
      dpotrf_( "L", &m, A, &m, &info );   
      
      get_ticks(ticksAFT); 
      ticks = ticksAFT - ticksB4;
      printf(" iter(%d): %.5g\n", i+1, (double) m*m*m / (double) (3*ticks) );
   }


      
   //mat_print( A, m, m, m, "L" );

   /*
   double one = 1.0, minus = -1.0;
   dgemm_( "N", "T", &m, &m, &m, &one, A, &m, A, &m, &minus, A0, &m ); 
   mat_print( A0, m, m, m, "diff" );
   */
   
   return 0;
}



void mat_print( double *A, int n, int m, int ldA, char *text)
{
   int i, j;
   FILE *output;
   output = stdout;

   fprintf(output, "%s = [ ...\n", text );
   for( i=0; i<n; i++ ) //row
      {
	 for( j=0; j<m; j++ )  //column
	    fprintf(output, "%.15g ", A[i+j*ldA] );
	 fprintf(output, "; ...\n" );
      }
   fprintf(output, "];\n" );

   return;
}