! ----------------------------------------------------------------------
! Testing LU decomposition by NR and LAPACK
! NR ludcmp & lubksb, LAPACK dgesv & dgetrf & dgetrs & dgesvx
! ----------------------------------------------------------------------
! Command-line compilation:
! ----------------------------------------------------------------------
! Linux
! gfortran -O3 mNR.for testLU-simple-x.f90 -llapack (default LAPACK library)
! ifort -heap-arrays mNR.for testLU-simple-x.f90 -mkl (Intel MKL library)
! nvfortran -O3 mNR.for testLU-simple-x.f90 -lblas (Nvidia HPC SDK library)
! ----------------------------------------------------------------------
! Windows
! gfortran -O3 mNR.for testLU-simple-x.f90 -llibblas (PGI library .. a problem with dgesvx)
! ifort -heap-arrays -o a mNR.for testLU-simple-x.f90 -Qmkl (Intel MKL library)
! (call setvars.bat or ifortvars.bat)
! gfortran -O3 testLU-simple-x.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(:),b_orig(:),x(:),work(:)
integer,allocatable,target :: indx(:),ipiv(:),iwork(:)
real(wp) :: d,r,c,rcond,ferr,berr
integer :: np,n,mp,m,nrhs,lda,ldb,ldx,info
character(1) :: trans,fact,equed
integer,parameter :: NR_GAUSSJ=1,NR_LU=2,NR_LUI=21,LAPACK=3,LAPACK2=4,LAPACKX=5
integer,parameter :: switch=LAPACKX

np=1000
n=np
mp=1
m=mp
allocate (a(np,np),b(np),x(np),indx(n),ipiv(n),work(n*4),iwork(n))
call random_number(a)
a_orig=a
b=1.
b_orig=b
select case (switch)
case (NR_GAUSSJ)
  call gaussj(a,n,np,b,m,mp)
  x=b
case (NR_LU,NR_LUI)
  call ludcmp(a,n,np,indx,d)
  call lubksb(a,n,np,indx,b)
  x=b
  if (switch==NR_LUI) call mprove(a_orig,a,n,np,indx,b_orig,x)
case (LAPACK)
  nrhs=m
  lda=np
  ldb=np
  call dgesv(n,nrhs,a,lda,ipiv,b,ldb,info)
  if (info/=0) print *,'dgesv info: ',info
  x=b
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)
  if (info/=0) print *,'dgetrf info: ',info
  call dgetrs(trans,n,nrhs,a,lda,ipiv,b,ldb,info)
  if (info/=0) print *,'dgetrs info: ',info
  x=b
case (LAPACKX)
  nrhs=m
  lda=np
  ldb=np
  ldx=np
  fact='N'  ! fact: 'N' no equilibration, 'E' equilibration
  trans='N' ! op(A)*X=B, trans: 'N' for no op, 'T' for transpose, 'C' for conjugate transpose
  equed='N' ! form of equilibration: 'N' when fact='N'
  call dgesvx(fact,trans,n,nrhs,a_orig,lda,a,lda,ipiv,equed,r,c,b,ldb,x,ldx,rcond,ferr,berr,work,iwork,info)
  print '(a,i0,2e20.12)','dgesvx info,rcond,ferr: ',info,rcond,ferr
case default
  print *,'Bad value of SWITCH.'
end select
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