//=================================================================C
//  Test Program of C=A*B for 4-Case                               C
//-----------------------------------------------------------------C
//    Written by Yasunori Ushiro ,  2007/05/30                     C
//        ( Tokyo Polytechnic University )                         C
//     後 保範（東京工芸大学）                                     C
//=================================================================C
#include <stdio.h> 
#include <math.h>
#include <time.h>
// Global Define
#define ND  1000
double A[ND][ND], B[ND][ND], C[ND][ND] ;
FILE   *FT1 ;
// Function Definition
void MULT1(int N, int M, int L) ;
void MULT2(int N, int M, int L) ;
void MULT3(double AA[][ND], double BB[][ND], double CC[][ND]
           , int N, int M, int L, int ID) ;
void MULT4(int N, int M, int L, int MB) ;
void SETAB(int N) ;
// Main Program
void main ()
{ int     i, j, k, N, M, L, MB, IER, T1, T2 ;
  double  CPU, FLOP, SUM, SUM1 ; 
  FT1 = fopen("C-MULT.txt","w") ;
  T1  = clock() ;                    // Start Time
//  Size Input
  printf("Type In  N(Matrix Size) and MB(Block Size) \n") ;
  scanf("%d %d",&N,&MB) ; 
  if(N > ND) { N = ND ; }
  if(MB > N) { MB = N ; }
  fprintf(FT1,"Matrix Multiplication (C=A*B)  N=%d  MB=%d \n",N,MB) ;
  fprintf(FT1,"Type IER  Time(s)  MFLOPS \n") ;
  printf("Type IER  Time(s)  MFLOPS \n") ; 
  M = N ;
  L = N ;
//  Loop for All Type of C=A*B
  for (k=1; k<=4; k++) {
//  Set A,B
    SETAB(N) ;
    T1 = clock() ;
    if(k <= 2) {
      if(k <= 1) {
//  JKI-Type Multiplication
        MULT1(N, M, L) ;  }
      else {
//  IKJ-Type Multiplication
        MULT2(N, M, L) ;  }
      } 
    else {
      if(k <= 3) {
//   Unrolling-Type Multiplication
        MULT3(A, B, C, N, M, L, 0) ;  }
      else {
//   Cache Blocking=type Multiplication
        MULT4(N, M, L, MB) ;  }
      }
    T2 = clock() ;
//  Output
    CPU = (double)(T2 - T1)/CLOCKS_PER_SEC ;
    FLOP = 0.0 ;
    if(CPU > 0.0) { FLOP = 2.0*N*N*N*1.0e-6/CPU ; }
    SUM = 0.0 ;
    for (i=0; i<N; i++) {
      for (j=0; j<N; j++) {
        SUM = SUM + C[i][j] ; }
      }
    if(k == 1) {
      SUM1 = SUM ;
      IER  = 0 ; }
    else {
      IER  = 1 ; 
      if(SUM == SUM1) { IER = 0 ; }
     }
    fprintf(FT1,"%3d %3d %8.1f %8.1f \n",k,IER,CPU,FLOP) ;   
    printf("%3d %3d %9.2f %8.1f \n",k,IER,CPU,FLOP) ;   
   }
 }
//=================================================================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 ;  }
    }
  }