! Linpack in Fortran 2003
!
! To run on F90/F95, remove the command line argument parsing part
!
! MJ Rutter, Nov 2015

program linpack

integer :: n,ic1,ic2,itick,info
integer, allocatable :: ipvt(:)
real(kind(1d0)), allocatable :: a(:,:),b(:),x(:),a2(:,:),b2(:),r(:)
real(kind(1d0)) :: eps, t1, ops,norma, resid, residn
character(len=20) :: arg

n=2000

! If using F95, start deleting here
if (command_argument_count().eq.1) then
  call get_command_argument(1,arg)
  read(arg,*) n
endif
! If using F95, end deleting here

eps=2.2e-16
ops=(2.0d0*n)*n*(n/3.0d0)+(2.0d0*n)*n

allocate(a(n,n),a2(n,n),b(n),b2(n),x(n),ipvt(n),r(n))

write(*,*)'Running Linpack with array size of ',n

call random_number(a)
a=a-0.5d0
a2=a

b=sum(a,dim=2)
b2=b

norma=maxval(abs(a))

call system_clock(ic1)

call dgesv(n,1,a,n,ipvt,b,n,info)

call system_clock(ic2,itick)

t1=real(ic2-ic1,kind(1d0))/itick

! restore over-written variables

a=a2
x=b
b=b2

! compute a residual to verify results

r=matmul(a,x)-b

resid=maxval(abs(r))
normx=maxval(abs(x))

residn=resid/(n*norma*normx*eps)

write(*,*)'Residual is ',resid
write(*,*)'Normalised residual is ',residn
write(*,*)'Machine epsilon is ',eps
write(*,*)'x(1)-1 is ',x(1)-1
write(*,*)'x(n)-1 is ',x(n)-1

write(*,*)'Time is ',t1
write(*,*)'MFLOPS: ',ops*1d-6/t1

end


