#include <stdio.h>
#include <math.h>
// Global define
#define ND  1000
extern double A[ND][ND], B[ND][ND], C[ND][ND] ;
// Function define
void MULT3(double *AA, double *BB, double *CC,
           int N, int M, int L, int ID) ;
int  min(int N1, int N2) ;
//=================================================================C
void MULT4(int N, int M, int L, int MB)
//=================================================================C
//  Real-Dense Matrix Multiplication                               C
//   C = A*B  with Cache Blocking                                  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
//    MB        I*4, In,  Blocking Size                            C
//-----------------------------------------------------------------C
//    Written by Yasunori Ushiro ,  2007/05/30                     C
//        ( Tokyo Polytechnic University )                         C
//     後 保範（東京工芸大学）                                     C
//=================================================================C
{ int   i, j, k, NS, MS, LS, ID ;   
//   C=A*B with Blocking
  for (i=0; i<N; i=i+MB) { 
    NS = min(MB, N-i) ;
    for (j=0; j<M; j=j+MB) { 
      MS = min(MB, M-j) ;
      ID = 0 ;
      for (k=0; k<L; k=k+MB) { 
        LS = min(MB, L-k) ;
        MULT3(&A[i][k], &B[k][j], &C[i][j], NS, MS, LS, ID) ;
        ID = 1 ;  }
    }  }
  }
//=================================================================C
int min(int N1, int N2)
{ int  MINV ;
  if(N1 <= N2) { MINV = N1 ; }
  else  { MINV = N2 ; }
  return (MINV) ;
 }