#include <stdio.h>
#include <math.h>
#define ND  1000
#define A(i,j)  A[i][j]
#define B(i,j)  B[i][j]
#define C(i,j)  C[i][j]
//=================================================================C
void MULT3(double A[][ND], double B[][ND], double C[][ND],
           int N, int M, int L, int ID)
//=================================================================C
//  Real-Dense Matrix Multiplication                               C
//   C = A*B  with i-k Unrolling in IKJ-Type                       C
//-----------------------------------------------------------------C
//    A[N][ND]  R*8, In,  First Input Matrix (A)                   C
//    B[L][ND]  R*8, In,  Second Input Matrix (B)                  C
//    C[N][ND]  R*8, Out, Output Matrix (C)                        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
//    ID        I*4, In,  ID=0,  Initial Set C=0                   C
//                          =1,  No Initial Set                    C
//-----------------------------------------------------------------C
//    Written by Yasunori Ushiro ,  2007/05/30                     C
//        ( Tokyo Polytechnic University )                         C
//     後 保範（東京工芸大学）                                     C
//=================================================================C
{ int     i, j, k, NE, LE, N2, L2 ;
  double  A1, A2, A3, A4 ;
//  Test of Even or Odd
  NE = N % 2 ;
  LE = L % 2 ;
  N2 = N - NE ;
  L2 = L - LE ;
//  Clear C
  if(ID == 0) {
    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<N2; i=i+2) {
    for (k=0; k<L2; k=k+2) {
      A1 = A(i,k) ;
      A2 = A(i,k+1) ;
      A3 = A(i+1,k) ;
      A4 = A(i+1,k+1) ;
      for (j=0; j<M; j++) {
        C(i,j)   = C(i,j)   + A1*B(k,j) + A2*B(k+1,j) ;
        C(i+1,j) = C(i+1,j) + A3*B(k,j) + A4*B(k+1,j) ;  }
    }  }
//  N is Odd
  if(NE == 1) {
    i = N - 1 ;
    for (k=0; k<L2; k=k+2) {
      A1 = A(i,k) ;
      A2 = A(i,k+1) ;
      for (j=0; j<M; j++) {
        C(i,j) = C(i,j) + A1*B(k,j) + A2*B(k+1,j) ;  }
    }  }
//  L is Odd
  if(LE == 1) {
    k = L - 1 ;
    for (i=0; i<N2; i=i+2) {
      A1 = A(i,k) ;
      A3 = A(i+1,k) ;
      for (j=0; j<M; j++) {
        C(i,j)   = C(i,j)   + A1*B(k,j) ;
        C(i+1,j) = C(i+1,j) + A3*B(k,j) ;  }
    }  } 
//  M and L are Odd
  if(NE == 1 & LE == 1) {
    i = N - 1 ;
    k = L - 1 ;
    A1 = A(i,k) ;
    for (j=0; j<M; j=j++) {
      C(i,j) = C(i,j) + A1*B(k,j) ;  }
   }
 }