m4_changequote({,})m4_define( {matmul},
{C
C----------------------------------------------------------------------
C
      SUBROUTINE MATMUL_$1$2$3(N, A, B, C)
      REAL*8 A(N,N), B(N,N), C(N,N)
      REAL WALLTIME, MFLOPS, ETIME, TARRAY(2), TIME0, TIME1, TIME2
*
      NTIMES = 100000000/(N*N*(2*N-1))
      IF ( NTIMES .EQ. 0 ) NTIMES = 1
*
      DO J = 1, N
      DO I = 1, N
         A(I,J) = 1.0
         B(I,J) = 1.0
      ENDDO
      ENDDO
*
      TIME0 = ETIME(TARRAY)
      DO L = 1, NTIMES
         DO J = 1, N
            DO I = 1, N
               C(I,J) = 0.0
            ENDDO
         ENDDO
         DO $1 = 1, N
            DO $2 = 1, N
               DO $3 = 1, N
                  C(I,J) = C(I,J) + A(I,K) * B(K,J)
               ENDDO
            ENDDO
         ENDDO
         CALL DUMMY(C)
      ENDDO
      TIME1 = ETIME(TARRAY)
      DO I = 1, NTIMES
         CALL DUMMY(C)
      ENDDO
      TIME2 = ETIME(TARRAY)
C
      WALLTIME = ((TIME1 - TIME0) - (TIME2 - TIME1)) / REAL(NTIMES)
      MFLOPS = (N*N*(2*N-1)) / (WALLTIME * 1.0E+06)
*
      PRINT *, '$1, $2, $3:  MFLOPS=', MFLOPS, ',   C(1,1)=', C(1,1)
      RETURN
      END
} )

matmul(I, J, K)
matmul(I, K, J)
matmul(J, I, K)
matmul(J, K, I)
matmul(K, I, J)
matmul(K, J, I)

c
c----------------------------------------------------------------------
c
      program matmul_time
      parameter (size=1025)
      real*8 a(size, size), b(size, size), c(size, size), arow(size)
c
      print *, 'Enter N: '
      read(*,*) n
      if ( n.gt.size ) then
         print *, 'N too large.'
         stop
      endif
*
      call matmul_ijk(n, a, b, c)
      call matmul_ikj(n, a, b, c)
      call matmul_jik(n, a, b, c)
      call matmul_jki(n, a, b, c)
      call matmul_kij(n, a, b, c)
      call matmul_kji(n, a, b, c)
*
      stop
      end
c
c----------------------------------------------------------------------
c
      subroutine dummy(c)
      double precision c
      return
      end