c This file includes modules for LU factoration and
c forward and backward substitution for banded matrices.

c nband: the bandwidth
c neqn:  the dimension of the matrix
c Note:  For single precision, replace "real*8" by "real*4".


c=======================================================================
      subroutine lufac(nband,neqn,level,A)
c=======================================================================
      implicit none
      integer nband,neqn,level
      real*8 A(-nband:nband,neqn)

      integer i,j,k,kb,kpi
      real*8  zero,pick,mult

      if(level.ge.1) write(*,'("LUFAC: level=",i1)') level
      if(level.ge.2) print*,"bandwidth=",nband," neqn=",neqn

      zero=0.0d0

       do 10 k=1,neqn
          kb = min(neqn-k,nband)
          do 20 i= 1, kb
             kpi = k+i
             pick= A(-i,kpi)
             if(pick.eq.zero) goto 20
             mult= pick/A(0,k)
             A(-i,kpi) = mult
             do j= 1, kb
                A(j-i,kpi)= A(j-i,kpi)-mult*A(j,k)
             end do
 20       continue
 10    continue

      return
      end

c=======================================================================
      subroutine substit(mx,neqn,level,A,x)
c=======================================================================
      implicit none
      integer mx,neqn,level
      real*8 A(-mx:mx,neqn),x(neqn)

      integer j,k,kb

      if(level.ge.1) write(*,'("SUBSTIT: level=",i1)') level
      if(level.ge.2) print*,"bandwidth=",mx," neqn=",neqn
 
      do k=2,neqn
         kb   = min(k-1,mx)
         do j= 1, kb
            x(k) = x(k) - A(-j,k)*x(k-j)
         end do
      end do

      x(neqn) = x(neqn)/A(0,neqn)
      do k= neqn-1, 1, -1
         kb= min(neqn-k,mx)
         do j= 1,kb
            x(k) = x(k) - A(j,k)*x(k+j)
         end do
         x(k) = x(k)/A(0,k)
      end do  

      return
      end