#include <stdio.h>
//================================================================C
// Test Program of Matrix-multiplication (C=A*B)                  C
//   with ikj-loop and (1x1,2x2)-unrolling                        C
//   by Single processor on Windows                               C
//----------------------------------------------------------------C
//   Written by Yasunori Ushiro,   2009/04/22                     C
//       ( Tokyu Kougei University )                              C
//================================================================C
#include <time.h>
#include <math.h>
#define n 5000  // Matrix size
#define lp 1    // Number of loops for use time
double A[n][n], B[n][n], C[n][n];
void set()
{ int i, j;
  double D;
  D = 1.0/n;
  for (i=0; i<n; i++)
    { for (j=0; j<n; j++)
        { A[i][j] = 1.0 + (2.0*i + j)*D;
          B[i][j] = 1.0 + (2.0*i - j)*D;
        }
    }
}
double sum()
{ double s;
  int i, j;
  s = 0.0;
  for (i=0; i<n; i++)
    { for (j=0; j<n; j++)
        { s += C[i][j]; }
    }
 return s; }
//-----------------------------------------------------C
// C=A*B, with ikj-loop and no-unrolling
int mult1()
{ int i, j, k;
  for (i=0; i<n; i++)
    { for (j=0; j<n; j++)
        { C[i][j] = 0.0; }
    }
  for (i=0; i<n; i++)
    { for (k=0; k<n; k++)
        { for (j=0; j<n; j++)
            { C[i][j] += A[i][k]*B[k][j]; }
        }
    }
return 0; }
//-----------------------------------------------------C
// C=A*B, with ikj-loop and 2x2-unrolling
int mult2()
{ int i, j, k;
  for (i=0; i<n; i++)
    { for (j=0; j<n; j++)
        { C[i][j] = 0.0; }
    }
  for (i=0; i<n; i+=2)
    { for (k=0; k<n; k+=2)
        { for (j=0; j<n; j++)
            { C[i][j]   += A[i][k]*B[k][j]
                         + A[i][k+1]*B[k+1][j];
              C[i+1][j] += A[i+1][k]*B[k][j]
                         + A[i+1][k+1]*B[k+1][j];
            }
        }
    }
return 0; }
//-----------------------------------------------------C
int main()
{ FILE *f1;
  int k, id, kj;
  long t1, t2;
  double use, Gf, gflops, K=1000.0, SG, s;
  SG = CLOCKS_PER_SEC;
  f1 = fopen("multpc.txt","w");
  Gf = 2.0*lp*(n/K)*(n/K)*(n/K);
  printf("Matrix multiplication, n=%d, lp=%d\n",n,lp);
  fprintf(f1,"Matrix multiplication, n=%d, lp=%d\n",n,lp);
  set();
// ikj-loop and no-unrolling 
  id = 0;
  t1 = clock();
  for (k=0; k<lp; k++)
     { id = mult1(); }
  t2 = clock();
  use = (double)(t2 - t1)/SG;
  gflops = Gf/use;
  s = sum();
  printf("no-unrolling,  %6.2f sec, %6.2f GFLOPS, %e\n",use,gflops,s);
  fprintf(f1,"no-unrolling,  %6.2f sec, %6.2f GFLOPS\n",use,gflops);
// ikj-loop and 2X2-unrolling 
  t1 = clock();
  for (k=0; k<lp; k++)
    { id = mult2(); }
  t2 = clock();
  use = (double)(t2 - t1)/SG;
  gflops = Gf/use;
  s = sum();
  printf("2x2-unrolling, %6.2f sec, %6.2f GFLOPS, %e\n",use,gflops,s);
  fprintf(f1,"2x2-unrolling, %6.2f sec, %6.2f GFLOPS\n",use,gflops);
  return 0; }