! ----------------------------------------------------------------------
! Testing LU decomposition of band diagonal matrices by NR and LAPACK
! NR bandec & banbks & banmul, LAPACK dgbtrf & dgbtrs & dgbmv
! ----------------------------------------------------------------------
! Command-line compilation:
! ----------------------------------------------------------------------
! Linux
! gfortran -O3 mNR.for testBand-simple.f90 -lblas -llapack (default LAPACK library)
! ifort -heap-arrays mNR.for testBand-simple.f90 -mkl (Intel MKL library)
! nvfortran -O3 mNR.for testBand-simple.f90 -lblas -llapack (Nvidia HPC SDK library)
! ----------------------------------------------------------------------
! Windows
! gfortran -O3 mNR.for testBand-simple.f90 -llibblas -lliblapack (PGI library ... a linking problem)
! ifort -heap-arrays -o a mNR.for testBand-simple.f90 -Qmkl (Intel MKL library)
! (call setvars.bat or ifortvars.bat)
! gfortran -O3 mNR.for testBand-simple.f90 "%MKLROOT%"\lib\intel64\mkl_rt.lib (Intel MKL library)
! ----------------------------------------------------------------------

program testBand
use mNR
implicit none
integer,parameter :: wp=8

real(wp),allocatable :: a(:,:),aorig(:,:),al(:,:),b(:),borig(:),x(:),work(:),ferr(:),berr(:)
integer,allocatable :: indx(:),iwork(:)
real(wp) d,anorm,rcond,dlange,alpha,beta
integer np,n,m1,m2,mp,mpl,nrhs,lda,ldb,ldab,info,i,incx,incy
character trans,norm
integer,parameter :: NR=1,LAPACK=2
integer,parameter :: switch=LAPACK

np=1000000
n=np
m1=60
m2=60
mp=m1+1+m2
mpl=m1
lda=mp
ldab=mpl+mp
ldb=np
nrhs=1

select case (switch)
case (NR)
  allocate (aorig(np,mp),a(np,mp),borig(np),b(np),al(np,mpl),x(np))
  allocate (indx(n))
case (LAPACK)
  allocate (aorig(lda,np),a(ldab,np),borig(ldb),b(ldb),x(ldb))
  allocate (indx(n))
case default
  print *,'Bad value of SWITCH.'
end select

call random_number(aorig)
borig=1.

do
select case (switch)
case (NR)
  a=aorig
case (LAPACK)
  a(m1+1:,:)=aorig
end select
b=borig

select case (switch)
case (NR)
  print *,'NR bandec & banbks'
  call bandec(a,n,m1,m2,np,mp,al,mpl,indx,d)
  call banbks(a,n,m1,m2,np,mp,al,mpl,indx,b)
  x=b
  print *,'NR banmul'
  call banmul(aorig,n,m1,m2,np,mp,x,b)
case (LAPACK)
  print *,'LAPACK dgbtrf & dgbtrs'
  trans='N'
  call dgbtrf(n,n,m1,m2,a,ldab,indx,info)
  if (info/=0) print *,'dgbtrf info =',info
  call dgbtrs(trans,n,m1,m2,nrhs,a,ldab,indx,b,ldb,info)
  if (info/=0) print *,'dgbtrs info =',info
  print *,'LAPACK dgbmv'
  x=b
  alpha=1.
  beta=0.
  incx=1
  incy=1
  call dgbmv(trans,n,n,m1,m2,alpha,aorig,lda,x,incx,beta,b,incy)
end select

print '(1x,a,10e20.12)','x =',x(1:3)
print '(1x,a,10e20.12)','b =',b(1:3)
enddo

end program