! ----------------------------------------------------------------------
! Testing LU decomposition of band diagonal matrices by LAPACK with NR interface
! ----------------------------------------------------------------------
! uncomment "use mBandec" below for NR solvers
! uncomment "use mBandecMKL" below for LAPACK solvers with NR interface
! ----------------------------------------------------------------------

module mConst
integer,parameter :: dp=8
  integer :: nmax=1000000,kmax=20,m1=60,m2=60
! integer :: nmax=3588,kmax=1,m1=6,m2=6
! integer :: nmax=3588,kmax=10000,m1=6,m2=6
end module

! ----------------------------------------------------------------------

module mBandecMKL
implicit none
integer,parameter,private :: dp=8
real(dp),allocatable,private :: ab(:,:)
integer,private :: ldab
integer,allocatable,private :: ipiv(:)

contains

SUBROUTINE bandec(a,n,m1,m2,np,mp,al,mpl,indx,d)
real(dp) :: a(:,:),al(:,:),d
integer n,m1,m2,np,mp,mpl,indx(:),m,info,i
intent(in) :: a,n,m1,m2,np,mp,al,mpl,indx
intent(out) :: d
m=m1+m2+1
ldab=m1+m
if (allocated(ab)) then
  if (size(ab,1)<ldab .or. size(ab,2)<n) deallocate (ab)
endif
if (.not.allocated(ab)) allocate (ab(ldab,n))
if (allocated(ipiv)) then
  if (size(ipiv)<n) deallocate (ipiv)
endif
if (.not.allocated(ipiv)) allocate (ipiv(n))
do i=1,m2+1
  ab(m1+i,m2+2-i:n)=a(1:n-m2-1+i,m+1-i)
enddo
do i=m2+2,m
  ab(m1+i,1:n+m2+1-i)=a(-m+m1+1+i:n,m+1-i)
enddo
ldab=size(ab,1)
call dgbtrf(n,n,m1,m2,ab,ldab,ipiv,info)
d=0.
END SUBROUTINE

SUBROUTINE banbks(a,n,m1,m2,np,mp,al,mpl,indx,b)
real(dp) :: a(:,:),al(:,:),b(n)
integer n,m1,m2,np,mp,mpl,indx(:),info
intent(in) :: a,n,m1,m2,np,mp,al,mpl,indx
intent(inout) :: b
integer,parameter :: nrhs=1
call dgbtrs('N',n,m1,m2,nrhs,ab,ldab,ipiv,b,n,info)
END SUBROUTINE

SUBROUTINE banmul(a,n,m1,m2,np,mp,x,b)
real(dp) :: a(np,mp),x(n),b(n)
integer n,m1,m2,np,mp,lda,i
intent(in) :: a,n,m1,m2,np,mp,x
intent(out) :: b
real(dp),allocatable :: aMKL(:,:)
real(dp),parameter :: alpha=1._8,beta=0._8
integer,parameter :: incx=1,incy=1
lda=m1+m2+1
allocate (aMKL(lda,n))
do i=1,m2+1
  aMKL(i,m2+2-i:n)=a(1:n-m2-1+i,lda+1-i)
enddo
do i=m2+2,lda
  aMKL(i,1:n+m2+1-i)=a(-lda+m1+1+i:n,lda+1-i)
enddo
call dgbmv('N',n,n,m1,m2,alpha,aMKL,lda,x,incx,beta,b,incy)
deallocate (aMKL)
END SUBROUTINE

SUBROUTINE evalRcond(a,n,m1,m2,rcond)
real(dp) :: a(:,:),anorm,rcond
integer n,m1,m2,info
intent(in) :: a,n,m1,m2
intent(out) :: rcond
character,parameter :: norm='1'  ! merge('1','I',.true.)
real(dp),allocatable :: work(:)
integer,allocatable :: iwork(:)
if (norm=='1') then
  anorm=maxval(sum(abs(a),1))
else
  anorm=maxval(sum(abs(a),2))
endif
allocate (work(3*n),iwork(n))
call dgbcon(norm,n,m1,m2,ab,ldab,ipiv,anorm,rcond,work,iwork,info)
deallocate (work,iwork)
END SUBROUTINE

end module

! ----------------------------------------------------------------------

module mBandec
implicit none

contains

SUBROUTINE bandec(a,n,m1,m2,np,mp,al,mpl,indx,d)
INTEGER m1,m2,mp,mpl,n,np,indx(n)
DOUBLE PRECISION d,a(np,mp),al(np,mpl),TINY
PARAMETER (TINY=1.d-20)
INTEGER i,j,k,l,mm
DOUBLE PRECISION dum
mm=m1+m2+1
if(mm.gt.mp.or.m1.gt.mpl.or.n.gt.np) print *,'bad args in bandec'
l=m1
do i=1,m1
  do j=m1+2-i,mm
    a(i,j-l)=a(i,j)
  enddo
  l=l-1
  do j=mm-l,mm
    a(i,j)=0.d0
  enddo
enddo
d=1.d0
l=m1
do k=1,n
  dum=a(k,1)
  i=k
  if(l.lt.n)l=l+1
  do j=k+1,l
    if(abs(a(j,1)).gt.abs(dum))then
      dum=a(j,1)
      i=j
    endif
  enddo
  indx(k)=i
  if(dum.eq.0.d0) a(k,1)=TINY
  if(i.ne.k)then
    d=-d
    do j=1,mm
      dum=a(k,j)
      a(k,j)=a(i,j)
      a(i,j)=dum
    enddo
  endif
  do i=k+1,l
    dum=a(i,1)/a(k,1)
    al(k,i-k)=dum
    do j=2,mm
      a(i,j-1)=a(i,j)-dum*a(k,j)
    enddo
    a(i,mm)=0.d0
  enddo
enddo
return
END SUBROUTINE

SUBROUTINE banbks(a,n,m1,m2,np,mp,al,mpl,indx,b)
INTEGER m1,m2,mp,mpl,n,np,indx(n)
DOUBLE PRECISION a(np,mp),al(np,mpl),b(n)
INTEGER i,k,l,mm
DOUBLE PRECISION dum
mm=m1+m2+1
if(mm.gt.mp.or.m1.gt.mpl.or.n.gt.np) print *,'bad args in banbks'
l=m1
do k=1,n
  i=indx(k)
  if(i.ne.k)then
    dum=b(k)
    b(k)=b(i)
    b(i)=dum
  endif
  if(l.lt.n)l=l+1
  do i=k+1,l
    b(i)=b(i)-al(k,i-k)*b(k)
  enddo
enddo
l=1
do i=n,1,-1
  dum=b(i)
  do k=2,l
    dum=dum-a(i,k)*b(k+i-1)
  enddo
  b(i)=dum/a(i,1)
  if(l.lt.mm) l=l+1
enddo
return
END SUBROUTINE

SUBROUTINE banmul(a,n,m1,m2,np,mp,x,b)
INTEGER m1,m2,mp,n,np
DOUBLE PRECISION a(np,mp),b(n),x(n)
INTEGER i,j,k
do i=1,n
  b(i)=0.d0
  k=i-m1-1
  do j=max(1,1-k),min(m1+m2+1,n-k)
    b(i)=b(i)+a(i,j)*x(j+k)
  enddo
enddo
return
END SUBROUTINE

SUBROUTINE evalRcond(a,n,m1,m2,rcond)
real(8) :: a(:,:),rcond
integer n,m1,m2
intent(in) :: a,n,m1,m2
intent(out) :: rcond
rcond=0.
END SUBROUTINE

end module

! ----------------------------------------------------------------------

program testBandecMKL
use mConst
  use mBandec
! use mBandecMKL
implicit none
real(dp),allocatable :: a(:,:),b(:),x(:),al(:,:),aa(:,:),bb(:)
integer,allocatable :: indx(:)
real(dp) :: d,rcond,t1,t2,tinit,tA,tB,tC
integer :: np,n,mp,mpl,k,input,id,i

np=nmax
n=np
mp=m1+1+m2
allocate (a(np,mp),b(np),x(np))
allocate (aa(np,mp),bb(np))
mpl=m1
allocate (al(np,mpl),indx(n))

tinit=0; tA=0; tB=0; tC=0

input=1
select case (input)
case (1)
  call random_number(aa)
case (2)
  call random_number(aa)
  aa(:,m1+1)=aa(:,m1+1)*nmax
case (3)
  open (newunit=id,file='matice')
  do i=1,nmax
  read (id,*) aa(i,:)
  enddo
  close (id)
end select
a=aa
bb=1.
b=bb

do k=1,kmax
  call cpu_time(t1)
  a=aa; b=bb
  call cpu_time(t2)
  tinit=tinit+(t2-t1)

  call cpu_time(t1)
  call bandec(a,n,m1,m2,np,mp,al,mpl,indx,d)
  call cpu_time(t2)
  tA=tA+(t2-t1)
  
  call cpu_time(t1)
  call banbks(a,n,m1,m2,np,mp,al,mpl,indx,b)
  x=b
  call cpu_time(t2)
  tB=tB+(t2-t1)
  
  call cpu_time(t1)
  call banmul(aa,n,m1,m2,np,mp,x,b)
  call cpu_time(t2)
  tC=tC+(t2-t1)
enddo

! call evalRcond(a,n,m1,m2,rcond)

print '(5e16.9)',x(1:4)
print '(5e16.9)',x(nmax-3:nmax)
print *,sum(abs(b-bb))/nmax
! print *,'rcond = ',rcond

deallocate (al,indx)
deallocate (aa,bb)

print '(5f16.9)',b(1:3)
print '(5f16.9)',b(nmax-2:nmax)
print '(4i8,4f8.3)',nmax,kmax,m1,m2,tinit,tA,tB,tC

end program