! -----------------------------------------------------------
! Example program
!
! Matrix multiplication using matmul vs. manual calc.
! Try out if you see any speed difference with statically allocated 
! tables and dynamically allocated ones in your system
! -----------------------------------------------------------
program test
  implicit none
  integer, parameter:: dbl=kind(1.d0),n=1000
  integer:: i,j,k
  real:: time0,time1
  real(dbl):: tmp
  real(dbl):: A(n,n),B(n,n),C(n,n)
!  real(dbl), allocatable,dimension(:,:):: A,B,C

!  allocate(A(n,n))
!  allocate(B(n,n))
!  allocate(C(n,n))


  print*,'Speed test of matrix multiplication: ',n,'x',n,' matrices'

  do i = 1, n
     do j = 1, n
        call random_number(A(i,j))
        call random_number(B(i,j))
     end do
  end do
  A = A*10.d0
  B = B*10.d0
  !
  ! matmul multiplication
  ! --------------------- 
  call cpu_time(time0)
  C=matmul(A,B)
  call cpu_time(time1)
  print*,'matmul timing ',time1-time0
  write(12,*) 'C(10,50) = ',C(10,50)
  !
  ! manual multiplication 
  ! --------------------- 
  call cpu_time(time0)

  do i = 1, n
     do j = 1, n
        tmp = 0.d0
        do k = 1, n
           tmp = tmp + A(i,k)*B(k,j)
        end do
        C(i,j) = tmp
     end do
  end do
  call cpu_time(time1)
  print*,'manual timing ',time1-time0
  write(12,*) 'C(10,50) = ',C(10,50)
  !
  ! BLAS multiplication
  ! -------------------
  call cpu_time(time0)
  call dgemm('N','N',n,n,n,1.d0,A,n,B,n,0.d0,C,n)
  call cpu_time(time1)
  print*,' BLAS timing ',time1-time0
  write(12,*) 'C(10,50) = ',C(10,50)

  end program test