program matmul_time
      parameter (size=1000)
      real*8 a(size, size), b(size, size), c(size, size), arow(size)
      real mflops0, mflops1
c
      print *, 'Enter N: '
      read(*,*) n
      if ( n.gt.size ) then
         print *, 'N too large.'
         stop
      endif
*
      call matmul0(n, a, b, c, mflops0)
      call matmul1(n, a, b, c, arow, mflops1)
*
      print *, 'N = ', N, '   MFLOPS0 = ', mflops0,
     +         '   MFLOPS1 = ', mflops1
      stop
      end
c
c----------------------------------------------------------------------
c
      subroutine matmul0(n, a, b, c, mflops)
      real*8 a(n,n), b(n,n), c(n,n)
      real walltime, mflops
      real etime
      real tarray(2), time0, time1, time2
c
      ntimes = 100000000/(n*n*(2*n-1))
      if ( ntimes .eq. 0 ) ntimes = 1
c
      do j = 1, n
      do i = 1, n
         a(i,j) = 1.0
         b(i,j) = 1.0
      enddo
      enddo
c
      time0 = etime(tarray)
      do l = 1, ntimes
         do j = 1, n
            do i = 1, n
               c(i,j) = 0.0
            enddo
         enddo
         do i = 1, n
            do j = 1, n
               do k = 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)
*
c     print *, "c(1,1) = ", c(1,1)
      return
      end
c
c----------------------------------------------------------------------
c
      subroutine matmul1(n, a, b, c, arow, mflops)
      real*8 a(n,n), b(n,n), c(n,n), arow(n)
      real walltime, mflops
      real etime
      real tarray(2), time0, time1, time2
c
      ntimes = 500000000/(n*n*(2*n-1))
      if ( ntimes .eq. 0 ) ntimes = 1
c
      do j = 1, n
      do i = 1, n
         a(i,j) = 1.0
         b(i,j) = 1.0
      enddo
      enddo
c
      time0 = etime(tarray)
      do l = 1, ntimes
         do j = 1, n
            do i = 1, n
               c(i,j) = 0.0
            enddo
         enddo
         do i = 1, n
            do k = 1, n
               arow(k) = a(i,k)
            enddo
            do j = 1, n
               do k = 1, n
                  c(i,j) = c(i,j) + arow(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)
c
c     print *, "c(1,1) = ", c(1,1)
      return
      end
c
c----------------------------------------------------------------------
c
      subroutine dummy(c)
      double precision c
      return
      end