subroutine ludcmpt(a, b, c, nmax)
      implicit double precision (o-z, a-h)
      dimension a(1), b(1), c(1)
      do n = nmax - 1, 1, -1
         if(dabs(b(n+1)) .lt. 1.d-30) then
            write(unit=*,fmt=*) n,' b = 0 in LUDCMPT '
            stop 111
         end if
         if(dabs(b(n+1)) .gt. 1.d+30) then
            write(unit=*,fmt=*) n+1,' b = inf in LUDCMPT '
            do m = nmax, 1, -1
               write(unit=*,fmt=*) m,a(m), b(m), c(m) 
            end do
            stop 112
         end if
         c(n) = c(n) / b(n + 1)
         b(n) = b(n) - c(n) * a(n + 1)
      end do
      return 
      end

      subroutine lubksbt(a, b, c, d, nmax)
      implicit double precision (o-z, a-h)
      dimension a(1), b(1), c(1), d(1), z(1000)
      z(nmax) = d(nmax)
      do n = nmax - 1, 1, -1
         z(n) = d(n) - (c(n) * z(n + 1))
      end do
      if(dabs(b(1)) .lt. 1.d-30) then
         write(unit=*,fmt=*) n,' b = 0 in LUBKSBT '
         stop 222 
      end if
      d(1) = z(1) / b(1)
      do n = 2, nmax, 1
         if(dabs(b(n)) .lt. 1.d-30) then
            write(unit=*,fmt=*) n,' b = 0 in LUBKSBT '
            stop 333 
         end if
         d(n) = (z(n) - (a(n) * d(n - 1))) / b(n)
      end do
      return 
      end