/*!
  \file   main.c
  \brief  Matrix multiplication implementation.

  \author Nikos Pitsianis
  \author Dimitris Floros
  \author Christos Choutouridis 8997 <cchoutou@ece.auth.gr>
  \date   2020-05-05
*/

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

#include "matmul.h"

extern double sqrt (double);

int cmpfunc (const void * a, const void * b) {
   double v =*(double*)a - *(double*)b;
   return (v < 0) ? -1 : (v > 0) ? 1 : 0;
}

double median (double* t, size_t n) {
   qsort ((void*)t, n, sizeof(t[0]), cmpfunc);
   return (n % 2) ? t[n/2] : (t[n/2] + t[n/2 -1]) /2;
}

double std_deviation (double* t, size_t n) {
   double av =0;
   for (size_t i=0 ; i<n ; ++i) {
      av += t[i];
   }
   av /= n;
   double s =0;
   for (size_t i=0 ; i<n ; ++i) {
      s += (t[i]-av)*(t[i]-av);
   }
   return sqrt (s/n);
}

mMult_ft multSelect (char* order) {
   if (! strcmp ((const char*)order, "ijk"))       return matrixMult_ijk;
   else if (! strcmp ((const char*)order, "ikj"))  return matrixMult_ikj;
   else if (! strcmp ((const char*)order, "jik"))  return matrixMult_jik;
   else if (! strcmp ((const char*)order, "jki"))  return matrixMult_jki;
   else if (! strcmp ((const char*)order, "kij"))  return matrixMult_kij;
   else if (! strcmp ((const char*)order, "kji"))  return matrixMult_kji;
   else                                            return matrixMult_ijk;
}

mMultBlock_ft multBlockSelect (char* order) {
   if (! strcmp ((const char*)order, "ijk"))       return matrixMult_ijk_block;
   else if (! strcmp ((const char*)order, "ikj"))  return matrixMult_ikj_block;
   else if (! strcmp ((const char*)order, "jik"))  return matrixMult_jik_block;
   else if (! strcmp ((const char*)order, "jki"))  return matrixMult_jki_block;
   else if (! strcmp ((const char*)order, "kij"))  return matrixMult_kij_block;
   else if (! strcmp ((const char*)order, "kji"))  return matrixMult_kji_block;
   else                                            return matrixMult_ijk_block;

}

/*!
 * A unit testing like main function to profile our code
 */
int main(int argc, char **argv) {
  struct timeval start, end;     /* time structs */
  double time[MAX_ITER] = {0.0}; /* execution time array in ms */
  float *A, *B, *C;              /* matrix declarations */
  mMult_ft mMult =NULL;
  mMultBlock_ft  mMultBlock =NULL;

  /* read matrix size (or use default) */
  if (argc < 3){
    fprintf( stderr,
          "Usage:\n"
          "%s n order [block], where \n"
          "   n: is the matrix size.\n"
          "   order: the loop order ex: ijk , jik, ...\n"
          "   block: Optional block size\n",
             argv[0]);
    exit(1);
  }
  int n = atoi( argv[1] );
  if (argc == 4)
      mMultBlock =  multBlockSelect(argv[2]);
  else
     mMult =  multSelect(argv[2]);

  /* initialize matrices */
  A = matrixInit( n );
  B = matrixInit( n );
  C = (float *) malloc( n*n*sizeof(float) );

   /* compute matrix multiplication */
   for (int it = 0; it < MAX_ITER; it++) {
      if (argc == 4) {
         gettimeofday(&start, NULL);
         mMultBlock( C, A, B, n, atoi(argv[3]) );
         gettimeofday(&end, NULL);
      }
      else {
         gettimeofday(&start, NULL);
         mMult( C, A, B, n );
         gettimeofday(&end, NULL);
      }
      time[it] =  (end.tv_sec  - start.tv_sec)  * 1000.0 +    /* sec to ms */
                  (end.tv_usec - start.tv_usec) / 1000.0;     /* us to ms */

      printf("Iter: %d Time: %f ms\n", it, time[it]);
   }

  /* we need to use the result -- verify it */
  for (int i = 0; i < n; i++) {             /* rows */
    for (int j = 0; j < n; j++) {           /* cols */
      float gold = 0;
      for (int k = 0; k < n; k++) {         /* accumulate products */
        gold += A[ sub2ind(i,k,n) ] * B[ sub2ind(k,j,n) ];
      }
      assert( (gold - C[sub2ind(i,j,n)]) < 1e-3 );
    }
  }

  // statistical data
  printf ("Median: %f [msec]\n", median (time, MAX_ITER));
  printf ("Std.Dev: %f [msec]\n", std_deviation(time, MAX_ITER));
}