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

    2007/05/30 東京工芸大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------

1. 概要

 実密行列を係数とする行列乗算C=A*Bを計算する。連続アクセス版の外側2個のループを
それぞれ、2でアンローリングしている。行列A,B,CのサイズはそれぞれNxL,LxM,NxMである。

2. プログラム

C=================================================================C
      SUBROUTINE MULT3(A,B,C,N,M,L,ND,ID)
C=================================================================C
C  Real-Dense Matrix Multiplication                               C
C   C = A*B  with j-k Unrolling in JKI-Type                       C
C-----------------------------------------------------------------C
C    A(ND,L)   R*8, In,  First Input Matrix (A) , Size N*L        C
C    B(ND,M)   R*8, In,  Second Input Matrix (B), Size L*M        C
C    C(ND,M)   R*8, Out, Output Matrix (C)      , Size N*M        C
C    N         I*4, In,  Matrix Size of A(Column), C(Column)      C
C    M         I*4, In,  Matrix Size of B(Row), C(Row)            C
C    L         I*4, In,  Matrix Size of A(Row), B(Column)         C
C    ND        I*4, In,  Array Size of A,B,C                      C
C    ID        I*4, In,  ID=0,  Initial Set C=0                   C
C                          =1,  No Initial Set                    C
C-----------------------------------------------------------------C
C    Written by Yasunori Ushiro ,  2007/05/30                     C
C        ( Tokyo Polytechnic University )                         C
C=================================================================C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(ND,L), B(ND,M), C(ND,M)
C  Test of Even or Odd
      ME = mod(M,2)
      LE = mod(L,2)
      M2 = M - ME
      L2 = L - LE
C  Clear C
      if(ID.eq.0) then
        do j=1,M
          do i=1,N
            C(i,j) = 0.0
          end do
        end do
      end if
C  C=C+A*B
      do j=1,M2,2
        do k=1,L2,2
          B1 = B(k,j)
          B2 = B(k+1,j)
          B3 = B(k,j+1)
          B4 = B(k+1,j+1)
          do i=1,N
            C(i,j)   = C(i,j)   + B1*A(i,k) + B2*A(i,k+1)
            C(i,j+1) = C(i,j+1) + B3*A(i,k) + B4*A(i,k+1)
          end do
        end do
      end do
C  M is Odd
      if(ME.eq.1) then
        j = M
        do k=1,L2,2
          B1 = B(k,j)
          B2 = B(k+1,j)
          do i=1,N
            C(i,j) = C(i,j) + B1*A(i,k) + B2*A(i,k+1)
          end do
        end do
      end if
C  L is Odd
      if(LE.eq.1) then
        k = L
        do j=1,M2,2
	    B1 = B(k,j)
          B3 = B(k,j+1)
          do i=1,N
            C(i,j)   = C(i,j)   + B1*A(i,k)
            C(i,j+1) = C(i,j+1) + B3*A(i,k)
          end do
        end do
      end if
C  M and L are Odd
      if(ME.eq.1 .or. LE.eq.1) then
	  j = M
        k = L
        B1 = B(k,j)
        do i=1,N
          C(i,j) = C(i,j) + B1*A(i,k)
        end do
      end if
C
      RETURN
      END