! ----------------------------------------------------------------------
! Testing LU decomposition by NR and LAPACK
! NR ludcmp & lubksb, LAPACK dgesv & dgetrf & dgetrs
! ----------------------------------------------------------------------
! Command-line compilation:
! ----------------------------------------------------------------------
! Linux
! gfortran -O3 mNR.for testLU-simple.f90 -llapack (default LAPACK library)
! ifort -heap-arrays mNR.for testLU-simple.f90 -mkl (Intel MKL library)
! nvfortran -O3 mNR.for testLU-simple.f90 -lblas (Nvidia HPC SDK library)
! ----------------------------------------------------------------------
! Windows
! gfortran -O3 mNR.for testLU-simple.f90 -llibblas (PGI library)
! ifort -heap-arrays -o a mNR.for testLU-simple.f90 -Qmkl (Intel MKL library)
! (call setvars.bat or ifortvars.bat)
! gfortran -O3 mNR.for testLU-simple.f90 "%MKLROOT%"\lib\intel64\mkl_rt.lib (Intel MKL library)
! ----------------------------------------------------------------------

program testLU
use mNR
implicit none
integer,parameter :: wp=8
real(wp),allocatable :: a(:,:),a_orig(:,:),b(:),x(:)
integer,allocatable :: indx(:),ipiv(:)
real(wp) :: d
integer :: np,n,mp,m,nrhs,lda,ldb,info
character(1) :: trans
integer,parameter :: NR_GAUSSJ=1,NR_LU=2,LAPACK=3,LAPACK2=4
integer,parameter :: switch=LAPACK

np=1000
n=np
mp=1
m=mp
allocate (a(np,np),b(np),indx(n),ipiv(n))
call random_number(a)
a_orig=a
b=1.
select case (switch)
case (NR_GAUSSJ)
  call gaussj(a,n,np,b,m,mp)
case (NR_LU)
  call ludcmp(a,n,np,indx,d)
  call lubksb(a,n,np,indx,b)
case (LAPACK)
  nrhs=m
  lda=np
  ldb=np
  call dgesv(n,nrhs,a,lda,ipiv,b,ldb,info)
  print *,'dgesv info: ',info
case (LAPACK2)
  nrhs=m
  lda=np
  ldb=np
  trans='N' ! op(A)*X=B, trans: 'N' for no op, 'T' for transpose, 'C' for conjugate transpose
  call dgetrf(n,n,a,lda,ipiv,info)
  print *,'dgetrf info: ',info
  call dgetrs(trans,n,nrhs,a,lda,ipiv,b,ldb,info)
  print *,'dgetrs info: ',info
case default
  print *,'Bad value of SWITCH.'
end select
x=b
print '(1x,a,10e20.12)','x =',x(1:3)
b=matmul(a_orig,x)  ! A.x=b
print '(1x,a,10e20.12)','b =',b(1:3)
end program