//=================================================================C
//  Test Program of C=A*B for 2-Case                               C
//-----------------------------------------------------------------C
//    Written by Yasunori Ushiro,   2007/04/08                     C
//        ( Kanagawa University )                                  C
//=================================================================C
#include <stdio.h> 
#include <math.h>
#include <time.h>
// Global Define
#define ND  2400
double A[ND][ND], B[ND][ND], C[ND][ND] ;
// Function Definition
//=================================================================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,   2003/08/28                     C
//        ( Hitachi Ltd. and 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 MULT1(int N, int M, int L)
//=================================================================C
//  Real-Dense Matrix Multiplication                               C
//   C = A*B  with 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,   2003/08/28                     C
//        ( Hitachi Ltd. and 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 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 ;  }
    }
  }
// Main Program
void main ()
{ int     i, j, k, N, M, L, IER, T1, T2, ID ;
  double  CPU, FLOP, SUM, SUM1 ; 
  T1  = clock() ;                    // Start Time
//  Size Input
LP1: printf("Type In  N (Matrix Size <= 2000) \n") ;
  scanf("%d",&N) ; 
  if(N > ND) { N = ND ; }
  printf("Type   Time(s)  MFLOPS \n") ; 
  M = N ;
  L = N ;
//  Loop for All Type of C=A*B
  for (k=1; k<=2; k++) {
//  Set A,B
    SETAB(N) ;
    T1 = clock() ;
    if(k <= 1) {
//  JKI-Type Multiplication
      MULT1(N, M, L) ;  }
    else {
//  IKJ-Type Multiplication
      MULT2(N, M, L) ;  }
    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 ; }
     }
    printf("%3d %9.2f %8.1f \n",k,CPU,FLOP) ;   
   }
  printf("Type in (1:continue, 0:stop) \n") ;
  scanf("%d",&ID) ;
  if(ID==1) goto LP1;
 }