! If using gfortran, force a call to a linked BLAS library from the matmaul()
! intrinsic by using the following flags:
!     -fexternal-blas -fblas-matmul-limit=1
! (See gfortran man page for details.)
!
! Otherwise replace matmul with a BLAS equivalent (e.g., dgemm).

implicit none

character(len=50) :: arg
integer :: n, iostat, i, j
real(8), dimension(:,:), allocatable :: a, b

n = 1000
if (command_argument_count() > 0) then
   call get_command_argument(1, arg)
   read(arg,*,iostat=iostat) n
   if (iostat /= 0) error stop 'First argument must be an integer n'
   if (n <= 0) error stop 'n must be > 0'
endif

allocate(a(n,n), b(n,n))

a = reshape([((i+j, j=1,n), i=1,n)], [n,n])
b = a

a = matmul(a, b)

end
