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,dclock
      integer i,j,ii,niter

      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 = dclock()
         do i = 1, niter
            call dgemm('n','n',m,n,k,alpha,a,ii,b,ii,beta,c,ii)
         enddo
         t2 = dclock()
         tt = t2 - t1
         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 = dclock()
         do i = 1, niter
            call dgemm('n','n',m,n,k,alpha,a,ii,b,ii,beta,c,ii)
         enddo
         t2 = dclock()
         tt = t2 - t1
         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 = dclock()
         do i = 1, niter
            call dgemm('n','n',m,n,k,alpha,a,ii,b,ii,beta,c,ii)
         enddo
         t2 = dclock()
         tt = t2 - t1
         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