#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/time.h>
#include <emmintrin.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;
            }
        }
      }
}


void  my_mm6_sse(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++) {
              __m128d aaa;
              __m128d aaa1;
              __m128d aaa2;
              __m128d aaa3;
              __m128d aaa4;
              __m128d xxx;
              __m128d www;

              ii = i + i1; jj = j + j1;
              tmp = *(c+ii*size + jj);
              aa = a + ii*size + k;
              ww = work + jj*size + k;

              aaa = _mm_loadl_pd(xxx, aa);
              aaa = _mm_loadh_pd(aaa, aa+1);

              aaa1 = _mm_loadl_pd(xxx, aa+2);
              aaa1 = _mm_loadh_pd(aaa1, aa+3);

              aaa2 = _mm_loadl_pd(xxx, aa+4);
              aaa2 = _mm_loadh_pd(aaa2, aa+5);

              aaa3 = _mm_loadl_pd(xxx, aa+6);
              aaa3 = _mm_loadh_pd(aaa3, aa+7);

              aaa4 = _mm_loadl_pd(xxx, aa+8);
              aaa4 = _mm_loadh_pd(aaa4, aa+9);


              www = _mm_loadl_pd(xxx, ww);
              www = _mm_loadh_pd(www, ww+1);
              aaa = _mm_mul_pd(aaa, www);

              www = _mm_loadl_pd(xxx, ww+2);
              www = _mm_loadh_pd(www, ww+3);
              aaa1 = _mm_mul_pd(aaa1, www);
              
              www = _mm_loadl_pd(xxx, ww+4);
              www = _mm_loadh_pd(www, ww+5);
              aaa2 = _mm_mul_pd(aaa2, www);

              www = _mm_loadl_pd(xxx, ww+6);
              www = _mm_loadh_pd(www, ww+7);
              aaa3 = _mm_mul_pd(aaa3, www);

              www = _mm_loadl_pd(xxx, ww+8);
              www = _mm_loadh_pd(www, ww+9);
              aaa4 = _mm_mul_pd(aaa4, www);

              aaa = _mm_add_pd(aaa, aaa4);
              aaa = _mm_add_pd(aaa, aaa3);
              aaa = _mm_add_pd(aaa, aaa2);
              aaa = _mm_add_pd(aaa, aaa1);

              aaa1 = _mm_loadl_pd(xxx, aa+10);
              aaa1 = _mm_loadh_pd(aaa1, aa+11);

              aaa2 = _mm_loadl_pd(xxx, aa+12);
              aaa2 = _mm_loadh_pd(aaa2, aa+13);

              aaa3 = _mm_loadl_pd(xxx, aa+14);
              aaa3 = _mm_loadh_pd(aaa3, aa+15);

              www = _mm_loadl_pd(xxx, ww+10);
              www = _mm_loadh_pd(www, ww+11);
              aaa1 = _mm_mul_pd(aaa1, www);

              www = _mm_loadl_pd(xxx, ww+12);
              www = _mm_loadh_pd(www, ww+13);
              aaa2 = _mm_mul_pd(aaa2, www);


              www = _mm_loadl_pd(xxx, ww+14);
              www = _mm_loadh_pd(www, ww+15);
              aaa3 = _mm_mul_pd(aaa3, www);

              aaa = _mm_add_pd(aaa, aaa1);
              aaa = _mm_add_pd(aaa, aaa2);
              aaa = _mm_add_pd(aaa, aaa3);
	      
              _mm_storel_pd(&t1, aaa);
              _mm_storeh_pd(&t0, aaa);

              tmp = t0+t1+tmp;
              *(c+ii*size + jj) = tmp;

            }
        }
      }
}

double ttt[2];

void  my_mm6_sse1(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++) {
              __m128d aaa;
              __m128d aaa1;
              __m128d aaa2;
              __m128d aaa3;
              __m128d aaa4;
              __m128d xxx;
              __m128d www;

              ii = i + i1; jj = j + j1;
              tmp = *(c+ii*size + jj);
              aa = a + ii*size + k;
              ww = work + jj*size + k;

              aaa = _mm_load_pd(aa);
              aaa1 = _mm_load_pd(aa+2);
              aaa2 = _mm_load_pd(aa+4);
              aaa3 = _mm_load_pd(aa+6);
              aaa4 = _mm_load_pd(aa+8);

              www = _mm_load_pd(ww);
              aaa = _mm_mul_pd(aaa, www);

              www = _mm_load_pd(ww+2);
              aaa1 = _mm_mul_pd(aaa1, www);
              
              www = _mm_load_pd(ww+4);
              aaa2 = _mm_mul_pd(aaa2, www);

              www = _mm_load_pd(ww+6);
              aaa3 = _mm_mul_pd(aaa3, www);

              www = _mm_load_pd(ww+8);
              aaa4 = _mm_mul_pd(aaa4, www);

              aaa = _mm_add_pd(aaa, aaa4);
              aaa = _mm_add_pd(aaa, aaa3);
              aaa = _mm_add_pd(aaa, aaa2);
              aaa = _mm_add_pd(aaa, aaa1);

              aaa1 = _mm_load_pd(aa+10);
              aaa2 = _mm_load_pd(aa+12);
              aaa3 = _mm_load_pd(aa+14);

              www = _mm_load_pd(ww+10);
              aaa1 = _mm_mul_pd(aaa1, www);

              www = _mm_load_pd(ww+12);
              aaa2 = _mm_mul_pd(aaa2, www);

              www = _mm_load_pd(ww+14);
              aaa3 = _mm_mul_pd(aaa3, www);

              aaa = _mm_add_pd(aaa, aaa1);
              aaa = _mm_add_pd(aaa, aaa2);
              aaa = _mm_add_pd(aaa, aaa3);
	      
              _mm_store_pd(ttt, aaa);

              tmp = ttt[0]+ttt[1]+tmp;
              *(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++) {
              __m128d aaa;
              __m128d aaa1;
              __m128d aaa2;
              __m128d aaa3;
              __m128d aaa4;
              __m128d www;

              ii = i + i1; jj = j + j1;
              tmp = *(c+ii*bk + jj);
              aa = a + ii*bk + k;
              ww = work + jj*bk + k;

	      /*
              _mm_loadl_pd(aaa, aa);
              _mm_loadh_pd(aaa, aa+1);

              _mm_loadl_pd(aaa1, aa+2);
              _mm_loadh_pd(aaa1, aa+3);

              _mm_loadl_pd(aaa2, aa+4);
              _mm_loadh_pd(aaa2, aa+5);

              _mm_loadl_pd(aaa3, aa+6);
              _mm_loadh_pd(aaa3, aa+7);

              _mm_loadl_pd(aaa4, aa+8);
              _mm_loadh_pd(aaa4, aa+9);


              _mm_loadl_pd(www, ww);
              _mm_loadh_pd(www, ww+1);
              aaa = _mm_mul_pd(aaa, www);

              _mm_loadl_pd(www, ww+2);
              _mm_loadh_pd(www, ww+3);
              aaa1 = _mm_mul_pd(aaa1, www);
              
              _mm_loadl_pd(www, ww+4);
              _mm_loadh_pd(www, ww+5);
              aaa2 = _mm_mul_pd(aaa2, www);

              _mm_loadl_pd(www, ww+6);
              _mm_loadh_pd(www, ww+7);
              aaa3 = _mm_mul_pd(aaa3, www);

              _mm_loadl_pd(www, ww+8);
              _mm_loadh_pd(www, ww+9);
              aaa4 = _mm_mul_pd(aaa4, www);

              aaa = _mm_add_pd(aaa, aaa4);
              aaa = _mm_add_pd(aaa, aaa3);
              aaa = _mm_add_pd(aaa, aaa2);
              aaa = _mm_add_pd(aaa, aaa1);


              _mm_loadl_pd(aaa1, aa+10);
              _mm_loadh_pd(aaa1, aa+11);

              _mm_loadl_pd(aaa2, aa+12);
              _mm_loadh_pd(aaa2, aa+13);

              _mm_loadl_pd(aaa3, aa+14);
              _mm_loadh_pd(aaa3, aa+15);

              _mm_loadl_pd(www, ww+10);
              _mm_loadh_pd(www, ww+11);
              aaa1 = _mm_mul_pd(aaa1, www);

              _mm_loadl_pd(www, ww+12);
              _mm_loadh_pd(www, ww+13);
              aaa2 = _mm_mul_pd(aaa2, www);


              _mm_loadl_pd(www, ww+14);
              _mm_loadh_pd(www, ww+15);
              aaa3 = _mm_mul_pd(aaa3, www);

              aaa = _mm_add_pd(aaa, aaa1);
              aaa = _mm_add_pd(aaa, aaa2);
              aaa = _mm_add_pd(aaa, aaa3);
	      */
              _mm_storel_pd(&tmp, aaa);
              _mm_storeh_pd(&t0, aaa);

              *(c+ii*bk + jj) = tmp+t0;
            }
        }
      }
}
 


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 = malloc(N*N*sizeof(double));
  B = malloc(N*N*sizeof(double));
  C = malloc(N*N*sizeof(double));
  W = 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 if (method == 8) 
      my_mm6_sse1(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);
   //   my_mm6_sse(N, B, A, C, W);
   for (i=0; i<N*N; i++)
     fprintf(fd, "%6.2lf\n", C[i]);
   fclose(fd);
 }
#endif
  
  return 0;
}