      program matmul_time
c
      integer i, j, k
      integer size, ntimes, m, n, p
      parameter (size=500)
      parameter (m=size,n=size,p=size)
      parameter (ntimes=2)
      real*8 a, b, c, arow
      dimension a(m,n), b(n,p), c(n,p), arow(n)

      integer l
      real*8 walltime, mflops
      integer hz, clock0, clock1, clock2
      integer t

      real*8 second
      real*8 time0, time1, time2

#if defined (OMPTIME)
      real*8 omp_get_wtime, omp_get_wtick
      print *, "OMP_Wtime resolution: ", omp_get_wtick()
#endif

#if defined (SYSCLK)
      call system_clock(count_rate=hz)
      print *, "system_clock resolution: ", real(1.d0/hz)
#endif
c
      do i = 1, m
      do j = 1, n
         a(i,j) = 1.0
      enddo
      enddo
      do i = 1, n
      do j = 1, p
         b(i,j) = 1.0
      enddo
      enddo

#if defined (SYSCLK)
! Elapsed time : Using F90 system_clock : Milli second
      call system_clock(count=clock0)

#elif defined (OMPTIME)
! Elapsed time : OpenMP build-in function: Micro second
      time0 = omp_get_wtime()
#else
! Elapsed time : SECOND() C function    : Micro second
      time0 = second()
#endif

      do l = 1, ntimes

#if defined MATMUL
! F90 Intrinsic function (MATMUL)
         c = matmul(a, b)

#elif defined BLAS
! AMD ACML(LAPACK) Library
         call dgemm('N','N', m, p, n, 1.0d0, a, m, b, n, 0.0d0, c, n)

#else
! Source coding
!$omp parallel default(shared)
!$omp do schedule(runtime)
         do j = 1, p
            do i = 1, m
               c(i,j) = 0.0
            enddo
         enddo
         do i = 1, m
!$omp do schedule(runtime)
            do ii = 1, n
               arow(ii) = a(i,ii)
            enddo
!$omp do schedule(runtime)
            do j = 1, p
               do k = 1, n
                  c(i,j) = c(i,j) + arow(k) * b(k,j)
               enddo
            enddo
         enddo
!$omp end parallel

#endif
         call dummy(c)
      enddo

#if defined (SYSCLK)
      call system_clock(count=clock1)
#elif defined (OMPTIME)
      time1 = omp_get_wtime()
#else
      time1 = second()
#endif

      do i = 1, ntimes
         call dummy(c)
      enddo

#if defined (SYSCLK)
      call system_clock(count=clock2)
#elif defined (OMPTIME)
      time2 = omp_get_wtime()
#else
      time2 = second()
#endif

#if defined (SYSCLK)
      t = (clock1 - clock0) - (clock2 - clock1)
      walltime = real(t) / real(hz) / real(ntimes)
      mflops = (m*p*(2*n-1)) / (walltime * 1.0d+06)
      print *,"Elapsed time =", walltime , " (sec)"
#else
      walltime = ((time1 - time0) - (time2 - time1)) / real(ntimes)
      mflops = (m*p*(2*n-1)) / (walltime * 1.0d+06)
      print *,"Elapsed time =", walltime ," (sec)"
#endif
c
      print *, "M =",M,", N =",N,", P =",P
      print *, "MFLOPS = ", mflops
      print *, "c(1,1) = ", c(1,1)
c
      end
c
      subroutine dummy(c)
      return
      end

