C=A*BのCプログラムNO.4

    2003/08/28 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------

1. 概要

 実密行列を係数とする行列乗算C=A*Bを計算する。計算機のキャッシュを有効に使用するため
ブロッキングしている。ブロックサイズは実行時に与える。NO.3のプログラムを使用している。
行列A,B,CのサイズはそれぞれNxL,LxM,NxMである。

2. プログラム

#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,   2003/08/28                     C
//        ( Hitachi Ltd. and Waseda University )                   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+1-i) ;
    for (j=0; j<M; j=j+MB) { 
      MS = min(MB, M+1-j) ;
      ID = 0 ;
      for (k=0; k<L; k=k+MB) { 
        LS = min(MB, L+1-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) ;
 }