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

  \author Nikos Pitsianis
  \author Dimitris Floros
  \date   2020-05-05
*/


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

#define MAX_ITER 10
#define sub2ind(i,j,n) (j) + (i)*(n)

/*
 * matrixMult - Matrix multiplication
 */
void matrixMult(float * const C,            /* output matrix */
                float const * const A,      /* first matrix */
                float const * const B,      /* second matrix */
                int const n) {              /* number of rows/cols */

  for (int i = 0; i < n; i++) {             /* rows */
    for (int j = 0; j < n; j++) {           /* cols */

      /* initialize output value */
      C[ sub2ind(i,j,n) ] = 0;

      for (int k = 0; k < n; k++) {         /* accumulate products */
        C[ sub2ind(i,j,n) ] +=
          A[ sub2ind(i,k,n) ] * B[ sub2ind(k,j,n) ];
      }
      
    }
  }
      
} // end function 'matrixMult'

/*
 * matrixInitAdd:
 *  Initialize matrix with random indices and return the matrix
 *  pointer.
 *                 
 */
float * matrixInit(int const n) {

  float *M = (float *) malloc( n*n*sizeof(float) );
  
  for (int i = 0; i < n; i++)              /* rows */
    for (int j = 0; j < n; j++)            /* cols */
      M[ sub2ind(i,j,n) ] = (float)rand()/(float)(RAND_MAX);

  return M;
  
} // end function 'matrixInitAdd'

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

  struct timeval start, end;    /* time structs */
  double time = 0.0;            /* execution time in ms */
  float *A, *B, *C;             /* matrix declarations */
  int n;                        /* matrix size */

  /* read matrix size (or use default) */
  if (argc != 2){
    fprintf( stderr, "Uasge:\n %s n\n   where n is the matrix size.\n",
             argv[0]);
    exit(1);
  }
  n = atoi( argv[1] );
  
  /* 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++) {
    gettimeofday(&start, NULL);
    matrixMult( C, A, B, n );
    gettimeofday(&end, NULL);

    time = ( (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);
  }

  /* we noeed 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 );
      
    }
  }
}