program matmul

  implicit none

  integer, parameter :: dp = selected_real_kind(14)
  integer, parameter :: n=1000 ! Matrix dimension
  real(kind=dp), dimension(n,n) :: A, B, C

  real(kind=dp) :: t0, t1 ! Timing variables

  integer :: i, j, k
  
  external dgemm

  ! Initialize matrices
  do j=1,n
     do i=1,n
        A(i,j) = i*j
        B(i,j) = 4*i - 6*j
     end do
  end do

  call cpu_time(t0)
  do j=1,n
     do i=1,n
        C(i,j) = 0.0
        do k=1,n
           C(i,j) = C(i,j) + A(i,k) * B(k,j)
        end do
     end do
  end do
  call cpu_time(t1)

  write(*,*) 'Elapsed CPU time with naive implementation', t1-t0

  call cpu_time(t0)
  call dgemm("n", "n", n, n, n, 1.0, A, n, B, n, 0.0, C, n)
  call cpu_time(t1)

  write(*,*) 'Elapsed CPU time with BLAS', t1-t0

end program matmul

