c        Copyright (C) 2009, The Portland Group, Incorporated.
c        Copyright (C) 2009, STMicroelectronics, Incorporated.
c        All rights reserved.
c
c          STMICROELECTRONICS, INCORPORATED PROPRIETARY INFORMATION
c   This software is supplied under the terms of a license agreement
c   or nondisclosure agreement with STMicroelectronics and may not be
c   copied or disclosed except in accordance with the terms of that
c   agreement.
c
c          This source code is intended only as a supplement to
c   STMicroelectronics Development Tools and/or on-line documentation.
c   See these sources for detailed information about
c   STMicroelectronics samples.
c
c          THIS CODE AND INFORMATION ARE PROVIDED "AS IS" WITHOUT
c   WARRANTY OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING BUT
c   NOT LIMITED TO THE IMPLIED WARRANTIES OF MERCHANTABILITY AND/OR
c   FITNESS FOR A PARTICULAR PURPOSE.
c

      program matmul
      real*8 a(2000,2000),b(2000,2000),c(2000,2000),alpha,beta
      integer m,n,k

      real*8 t1,t2,tt,mflops
      integer i,j,ii,niter

* Elapsed time function initialized -- SofTek
      integer ic, icr, icm
      real*8   t_ac
      real*8   second
      external second
      call system_clock(ic,icr,icm)
      t_ac= 1.0d0/real(icr,8)
      print '(1x,a,e10.5)','Time measurement accuracy : ',t_ac

      do j = 1, 2000
         do i = 1, 2000
            a(i,j) = 1.0d0
            b(i,j) = 1.0d0
            c(i,j) = 0.0d0
         enddo
      enddo

      alpha = 1.0d0
      beta = 0.0d0

      do ii = 10, 100, 10

         if (ii .le. 30) then
            niter = 10000
         else
            niter = 1000
         endif

         m = ii
         n = ii
         k = ii

         t1 = second()
         do i = 1, niter
            call dgemm('n','n',m,n,k,alpha,a,ii,b,ii,beta,c,ii)
         enddo
         t2 = second()
         tt = (t2 - t1)*t_ac
         mflops = 2.0d0 * dble(m) * dble(n) * dble(k) * dble(niter)
         mflops = mflops / 1000000.0d0
         mflops = mflops / tt

c..         print *, "Total Time(n): ",n,tt
c..         print *, "Mflops(n): ",n,mflops
         print *, n,mflops

      enddo


      do ii = 100, 1000, 100

         if (ii .le. 300) then
            niter = 100
         else
            niter = 10
         endif

         m = ii
         n = ii
         k = ii

         t1 = second()
         do i = 1, niter
            call dgemm('n','n',m,n,k,alpha,a,ii,b,ii,beta,c,ii)
         enddo
         t2 = second()
         tt = (t2 - t1)*t_ac
         mflops = 2.0d0 * dble(m) * dble(n) * dble(k) * dble(niter)
         mflops = mflops / 1000000.0d0
         mflops = mflops / tt

c..         print *, "Total Time(n): ",n,tt
c..         print *, "Mflops(n): ",n,mflops
         print *, n,mflops

      enddo

      niter = 1
      do ii = 1000, 2000, 100

         m = ii
         n = ii
         k = ii

         t1 = second()
         do i = 1, niter
            call dgemm('n','n',m,n,k,alpha,a,ii,b,ii,beta,c,ii)
         enddo
         t2 = second()
         tt = (t2 - t1)*t_ac
         mflops = 2.0d0 * dble(m) * dble(n) * dble(k) * dble(niter)
         mflops = mflops / 1000000.0d0
         mflops = mflops / tt

c..         print *, "Total Time(n): ",n,tt
c..         print *, "Mflops(n): ",n,mflops
         print *, n,mflops

      enddo

      end

      function second() result(rtime)
      implicit none
!
      integer :: ic,ir,im
      real(8) :: rtime
!
      call system_clock(ic,ir,im)
      rtime= real(ic,8)
!
      return
      end