!      ___________________________________________________
!      Single/Multiple Node Fortran 90 matrix-matrix multiplication
!      Donev Aleksandar, 11/13/99
!      ___________________________________________________

       PROGRAM MMul

       USE blas_kinds
       USE blas_keywords
       USE interface_gemm, ONLY: gemm

         INTEGER :: i,j,k
         INTEGER, PARAMETER :: n=500
         REAL(dp), DIMENSION(n,n) :: a,b,c
         REAL :: t(2),t0,t1

         DO i=1,n
           DO j=1,n
             a(i,j)=1.0_dp*i*j
             b(i,j)=1.0_dp*(i+j)
           END DO
         END DO
         c=0.0_dp  
    
         t0=ETIME(t)
         DO j=1,n
           DO i=1,n
             DO k=1,n
               c(i,j)=c(i,j)+a(i,k)*b(k,j)
             END DO
           END DO
         END DO

         t1=ETIME(t)
         WRITE(*,*) "C(34,67)=",c(34,67)
         WRITE(*,*) "Elapsed time for nested DO's =",t1-t0
         
         t0=ETIME(t)
         c=MATMUL(a,b)
         t1=ETIME(t)
         WRITE(*,*) "C(34,67)=",c(34,67)
         WRITE(*,*) "Elapsed time for intrinsic MMUL=",t1-t0

         t0=ETIME(t)
         call gemm(a,b,c)
         t1=ETIME(t)
         WRITE(*,*) "C(34,67)=",c(34,67)
         WRITE(*,*) "Elapsed time for BLAS gbmv=",t1-t0

       END PROGRAM MMul