//================================================================C
//  Speed Test Program of C=A*B for 2-Type                        C
//    MULT1: Serial access (IKJ Type) ; High speed type           C
//    MULT2: Stride access (JKI Type) ; Low speed type            C
//----------------------------------------------------------------C
//    Written by Yasunori Ushiro ,  2011/04/05                    C
//        ( JAMSTEC ES Adviser & Waseda University )              C
//================================================================C
#include <stdio.h> 
#include <math.h>
#include <time.h>
// Global Define
#define ND  2048
double A[ND][ND], B[ND][ND], C[ND][ND];
//================================================================C
void MULT1(int N, int M, int L)
//================================================================C
//  Real-Dense Matrix Multiplication                              C
//   C = A*B  with sirial access (IKJ-Type)                       C
//----------------------------------------------------------------C
//    A,B,C     Global define Matrix                              C
//    N         I*4, In,  Matrix Size of A(Column), C(Column)     C
//    M         I*4, In,  Matrix Size of B(Row), C(Row)           C
//    L         I*4, In,  Matrix Size of A(Row), B(Column)        C
//----------------------------------------------------------------C
//    Written by Yasunori Ushiro ,  2011/04/05                    C
//        ( JAMSTEC ES Adviser & Waseda University )              C
//================================================================C
{ int     i, j, k ;
  double  AW ;
//   Clear C
  for (i=0; i<N; i++) {
    for (j=0; j<M; j++) {
      C[i][j] = 0.0 ; }
   }
//   C=C+A*B
  for (i=0; i<N; i++) {
    for (k=0; k<L; k++) {
      AW = A[i][k] ;
      for (j=0; j<M; j++) {
        C[i][j] = C[i][j]  + AW*B[k][j] ; }
   }  }
 }
//================================================================C
void MULT2(int N, int M, int L)
//================================================================C
//  Real-Dense Matrix Multiplication                              C
//   C = A*B  with JKI-Type                                       C
//----------------------------------------------------------------C
//    A,B,C     Global define Matrix                              C
//    N         I*4, In,  Matrix Size of A(Column), C(Column)     C
//    M         I*4, In,  Matrix Size of B(Row), C(Row)           C
//    L         I*4, In,  Matrix Size of A(Row), B(Column)        C
//----------------------------------------------------------------C
//    Written by Yasunori Ushiro ,  2011/04/05                    C
//        ( JAMSTEC ES Adviser & Waseda University )              C
//================================================================C
{ int     i, j, k ;
  double  BW ;
//   Clear C
  for (j=0; j<M; j++) {
    for (i=0; i<N; i++) {
      C[i][j] = 0.0 ; }
   }
//   C=C+A*B
  for (j=0; j<M; j++) {
    for (k=0; k<L; k++) {
      BW = B[k][j] ;
      for (i=0; i<N; i++) {
        C[i][j] = C[i][j] + A[i][k]*BW ; }
   }  }
 }
//================================================================C
void SETAB(int N)
//----------------------------------------------------------------C
{ int  i, j ;
//  Set A,B
  for (i=0; i<N; i++) {
    for (j=0; j<N; j++) {
      A[i][j] = i + j + 2 ;
      B[i][j] = 2*i - j ;  }
    }
  }
//================================================================C
//  Speed Test Main Program of C=A*B for 2-Type                   C
//    MULT1: Serial access (IKJ Type) ; High speed type           C
//    MULT2: Stride access (JKI Type) ; Low speed type            C
//----------------------------------------------------------------C
//    Written by Yasunori Ushiro ,  2011/04/05                    C
//        ( JAMSTEC ES Adviser & Waseda University )              C
//================================================================C
// Main Program
int main ()
 { int N, M, L;
   long long int T1, T2;
   double  CPU, FLOP, SUM, SUM1 ; 
//  Size Input
   printf("Type In  N(Matrix Size) <= 2000\n") ;
   scanf("%d",&N) ; 
   if(N > ND) { N = ND ; }
   printf("Matrix Multiplication (C=A*B)  size=%d\n",N) ;
   M = N ;
   L = N ;
//  Loop for two Type of C=A*B
//  Set A,B
   SETAB(N);
//  IKJ-Type(Serial access) Multiplication
   T1 = clock();
   MULT1(N, M, L);
   T2 = clock();
   CPU = (double)(T2 - T1)/CLOCKS_PER_SEC;
   FLOP = 0.0;
   if(CPU > 0.0) { FLOP = 2.0*N*N*N*1.0e-6/CPU; }
   printf("Serial access; CPU time(s)=%7.2f, MFLOPS=%7.1f \n",CPU,FLOP); 
//  JKI-Type(Stride access) Multiplication
   T1 = clock();  
   MULT2(N, M, L); 
   T2 = clock();
   CPU = (double)(T2 - T1)/CLOCKS_PER_SEC;
   FLOP = 0.0;
   if(CPU > 0.0) { FLOP = 2.0*N*N*N*1.0e-6/CPU; }
   printf("Stride access; CPU time(s)=%7.2f, MFLOPS=%6.1f \n",CPU,FLOP); 
   return 0;
 }