/*
 * example program for CDA5125
 *
 * To compile: gcc -O3 my_mm.c
 *
 * To run: ./a.out matrix_size iter method
 * Example: to perform matrix multiplication of 1024x1024 matrices
 *          2 times using method 0 -- './a.out 1024 2 0'
 *           
 * Note: This program assumes that the matrix size is a multiple of 16
 *       and does not deal with the cases when the size is not a multiple
 *       of 16.
 */

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

struct timeval t1, t2;

void  naive_mm(int size, double *a, double *b, double *c) 
{
  int i, j, k;
    for (i=0; i<size; i++)
      for (j=0; j < size; j++)
        for (k=0; k<size; k++) {
          *(c+i*size + j) = *(c+i*size + j) + *(a+i*size + k) *  *(b +k*size + j);
	}
}

void  my_mm1(int size, double *a, double *b, double *c, double *work) 
{
  int i, j, k;

  for (i=0; i<size; i++)
    for (j=0; j<size; j++) *(work+i*size+j) = *(b + i+j*size);

    for (i=0; i<size; i++)
      for (j=0; j < size; j++)
        for (k=0; k<size; k++) {
          *(c+i*size + j) = *(c+i*size + j) + *(a+i*size + k) *  *(work +k+ size *j);
        }
}

void  my_mm2(int size, double *a, double *b, double *c, double *work) 
{
  int i, j, k;
  double tmp;
  for (i=0; i<size; i++)
    for (j=0; j<size; j++) *(work+i*size+j) = *(b + i+j*size);

    for (i=0; i<size; i++)
      for (j=0; j < size; j++) {
        tmp = *(c+i*size + j);
        for (k=0; k<size; k++) {
          tmp = tmp + *(a+i*size + k) *  *(work +k+ size *j);
        }
        *(c+i*size + j) = tmp;
      }
}


void  my_mm3(int size, double *a, double *b, double *c, double *work) 
{
  int i, j, k, i1, j1, k1, bk, ii, jj;
  double tmp;

  bk = 16;

  for (i=0; i<size; i++)
    for (j=0; j<size; j++) *(work+i*size+j) = *(b + i+j*size);

  for (i=0; i<size; i+=bk)
    for (j=0; j < size; j+=bk) {

      for (i1 = 0; i1<bk; i1++)
	for (j1 = 0; j1<bk; j1++) {
	  ii = i + i1; jj = j + j1;

	  tmp = *(c+ii*size + jj);
	  for (k=0; k<size; k++) {
	    tmp = tmp + *(a+ii*size + k) *  *(work +k+ size *jj);
	  }
	  *(c+ii*size + jj) = tmp;
	}
    }
}


void  my_mm4(int size, double *a, double *b, double *c, double *work) 
{
  int i, j, k, i1, j1, k1, bk, ii, jj, kk;
  double tmp;

  bk = 16;

  for (i=0; i<size; i++)
    for (j=0; j<size; j++) *(work+i*size+j) = *(b + i+j*size);

  //  for (i=0; i<size; i++)
  //  for (j=0; j<size; j++) *(work+i*size+j) = *(b + i+j*size);

    for (i=0; i<size; i+=bk)
      for (j=0; j < size; j+=bk) {
        for (k=0; k<size; k+=bk) {

          for (i1 = 0; i1<bk; i1++)
            for (j1 = 0; j1<bk; j1++) {
              ii = i + i1; jj = j + j1;
              tmp = *(c+ii*size + jj);
              for (k1=0; k1<bk; k1++) {
                kk = k + k1;
                tmp = tmp + *(a+ii*size + kk) *  *(work +kk+ size *jj);
              }
            
              *(c+ii*size + jj) = tmp;
            }
        }
      }
}


void  my_mm5(int size, double *a, double *b, double *c, double *work) 
{
  int i, j, k, i1, j1, k1, bk, ii, jj, kk;
  double tmp,t1, t2,t3,t4,t5,t6,t7,t8,t9, t10, t11, t12, t13, t14, t15, t0;
  double *aa, *ww;

  bk = 16;

  for (i=0; i<size; i++)
    for (j=0; j<size; j++) *(work+i*size+j) = *(b + i+j*size);

    for (i=0; i<size; i+=bk)
      for (j=0; j < size; j+=bk) {
        for (k=0; k<size; k+=bk) {

          for (i1 = 0; i1<bk; i1++)
            for (j1 = 0; j1<bk; j1++) {
              ii = i + i1; jj = j + j1;
              tmp = *(c+ii*size + jj);
              aa = a + ii*size + k;
              ww = work + jj*size + k;
            
              tmp = tmp + *(aa ) *  *(ww );
              tmp = tmp + *(aa +1) *  *(ww +1);
              tmp = tmp + *(aa +2) *  *(ww +2);
              tmp = tmp + *(aa +3) *  *(ww +3);
              tmp = tmp + *(aa +4) *  *(ww +4);
              tmp = tmp + *(aa +5) *  *(ww +5);
              tmp = tmp + *(aa +6) *  *(ww +6);
              tmp = tmp + *(aa +7) *  *(ww +7);
              tmp = tmp + *(aa +8) *  *(ww +8);
              tmp = tmp + *(aa +9) *  *(ww +9);
              tmp = tmp + *(aa +10) *  *(ww +10);
              tmp = tmp + *(aa +11) *  *(ww +11);
              tmp = tmp + *(aa +12) *  *(ww +12);
              tmp = tmp + *(aa +13) *  *(ww +13);
              tmp = tmp + *(aa +14) *  *(ww +14);
              tmp = tmp + *(aa +15) *  *(ww +15);
            
              *(c+ii*size + jj) = tmp;
            }
        }
      }
}



void  my_mm6(int size, double *a, double *b, double *c, double *work) 
{
  int i, j, k, i1, j1, k1, bk, ii, jj, kk;
  double tmp,t1, t2,t3,t4,t5,t6,t7,t8,t9, t10, t11, t12, t13, t14, t15, t0;
  double *aa, *ww;

  bk = 16;

  for (i=0; i<size; i++)
    for (j=0; j<size; j++) *(work+i*size+j) = *(b + i+j*size);

    for (i=0; i<size; i+=bk)
      for (j=0; j < size; j+=bk) {
        for (k=0; k<size; k+=bk) {

          for (i1 = 0; i1<bk; i1++)
            for (j1 = 0; j1<bk; j1++) {
              ii = i + i1; jj = j + j1;
              tmp = *(c+ii*size + jj);
              aa = a + ii*size + k;
              ww = work + jj*size + k;
            
              t0 =  *(aa  ) *  *(ww  );
              t1 = *(aa  +1) *  *(ww  +1);
              t2 = *(aa  +2) *  *(ww  +2);
              t3 = *(aa  +3) *  *(ww  +3);
              t4 =  *(aa  +4) *  *(ww  +4);
              t5 =  *(aa  +5) *  *(ww  +5);
              t6 =  *(aa  +6) *  *(ww  +6);
              t7 = *(aa  +7) *  *(ww  +7);
              t8 =  *(aa  +8) *  *(ww  +8);
              t9 = *(aa  +9) *  *(ww  +9);
              t10 = *(aa  +10) *  *(ww  +10);
              t11 = *(aa  +11) *  *(ww  +11);
              t12 = *(aa  +12) *  *(ww  +12);
              t13 = *(aa  +13) *  *(ww  +13);
              t14 = *(aa  +14) *  *(ww  +14);
              t15 = *(aa  +15) *  *(ww  +15);
            
              tmp = tmp + t0 + t1 + t2 + t3 + t4 + t5 + t6+t7
                +t8 + t9 + t10 + t11 + t12 + t13 + t14 + t15;

              *(c+ii*size + jj) = tmp;
            }
        }
      }
}


// this routine is incorrect yet.
// but it should have the same performance as block 
// order allocation.

void  my_mm7(int size, double *a, double *b, double *c, double *work) 
{
  int i, j, k, i1, j1, k1, bk, ii, jj, kk;
  double tmp,t1, t2,t3,t4,t5,t6,t7,t8,t9, t10, t11, t12, t13, t14, t15, t0;
  double *aa, *ww;

  bk = 16;

  for (i=0; i<size; i++)
    for (j=0; j<size; j++) *(work+i*size+j) = *(b + i+j*size);

  for (i=0; i<size; i++)
    for (j=0; j<size; j++) *(work+i*size+j) = *(b + i+j*size);

    for (i=0; i<size; i+=bk)
      for (j=0; j < size; j+=bk) {
        for (k=0; k<size; k+=bk) {

          for (i1 = 0; i1<bk; i1++)
            for (j1 = 0; j1<bk; j1++) {
              ii = i + i1; jj = j + j1;
              tmp = *(c+ii*bk + jj);
              aa = a + ii*bk + k;
              ww = work + jj*bk + k;
            
              t0 =  *(aa  ) *  *(ww  );
              t1 = *(aa  +1) *  *(ww  +1);
              t2 = *(aa  +2) *  *(ww  +2);
              t3 = *(aa  +3) *  *(ww  +3);
              t4 =  *(aa  +4) *  *(ww  +4);
              t5 =  *(aa  +5) *  *(ww  +5);
              t6 =  *(aa  +6) *  *(ww  +6);
              t7 = *(aa  +7) *  *(ww  +7);
              t8 =  *(aa  +8) *  *(ww  +8);
              t9 = *(aa  +9) *  *(ww  +9);
              t10 = *(aa  +10) *  *(ww  +10);
              t11 = *(aa  +11) *  *(ww  +11);
              t12 = *(aa  +12) *  *(ww  +12);
              t13 = *(aa  +13) *  *(ww  +13);
              t14 = *(aa  +14) *  *(ww  +14);
              t15 = *(aa  +15) *  *(ww  +15);
            
              tmp = tmp + t0 + t1 + t2 + t3 + t4 + t5 + t6+t7
                +t8 + t9 + t10 + t11 + t12 + t13 + t14 + t15;
              *(c+ii*bk + jj) = tmp;
            }
        }
      }
}
 


int
main( int argc, char *argv[])
{
  double *A, *B, *C, *W, *Z, *WORK;
  int  N, M, i;
  int method;

  if (argc < 3) {
    printf("Usage: a.out matrix_size iter opt_method\n");
    exit(0);
  }

  N= atoi(argv[1]);
  M = atoi(argv[2]);
  method = 0;
  if (argc >= 4) 
    method = atoi(argv[3]);

  srand48(100);

  A = (double *)malloc(N*N*sizeof(double));
  B = (double *)malloc(N*N*sizeof(double));
  C = (double *)malloc(N*N*sizeof(double));
  W = (double *)malloc(N*N*sizeof(double));
  for (i=0; i<N*N; ++i) {
    A[i] = drand48();
    B[i] = drand48();
  }

  
  gettimeofday(&t1, NULL);
  for (i=0; i<M; i++) {
    if (method == 0) 
      naive_mm(N,B,A,C);
    else if (method == 1) 
      my_mm1(N, B, A, C, W);
    else if (method == 2) 
      my_mm2(N, B, A, C, W);
    else if (method == 3) {
      printf("Don't want to show the results for this method.\n");
      exit(0);
    }
    //      my_mm3(N, B, A, C, W);
    else if (method == 4) 
      my_mm4(N, B, A, C, W);
    else if (method == 5) 
      my_mm5(N, B, A, C, W);
    else if (method == 6) 
      my_mm6(N, B, A, C, W);
    else if (method == 7) 
      my_mm7(N, B, A, C, W);
    else {
      printf("Method not supported.\n");
      exit(0);
    }
  }
  gettimeofday(&t2, NULL);

  printf("Time for the matrix multiplication using method %d is %d milliseconds\n",
	 method,
         (t2.tv_sec - t1.tv_sec)*1000 + 
         (t2.tv_usec - t1.tv_usec) / 1000);

#ifdef CHECK
 {
   FILE *fd;
   if ((fd = fopen("tmp333", "w")) == NULL) {
     printf("Cannot open tmp333\n"); exit(0);
   }

   naive_mm(N, B, A, C);   
   //   my_mm6(N, B, A, C, W);
   for (i=0; i<N*N; i++)
     fprintf(fd, "%6.2lf\n", C[i]);
   fclose(fd);
 }
#endif
  
  return 0;
}