C=A*BのCプログラムNO.3
2003/08/28 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
実密行列を係数とする行列乗算C=A*Bを計算する。連続アクセス版の外側2個のループを
それぞれ、2でアンローリングしている。行列A,B,CのサイズはそれぞれNxL,LxM,NxMである。
2. プログラム
#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, 2003/08/28 C
// ( Hitachi Ltd. and Waseda University ) 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) ; }
}
}