subroutine ludcmp(a, n, np, indx, d, ierr)
      implicit double precision (o-z, a-h)
      parameter (nmax = 1106, tiny = 1.0d-30, zero = 0.d0)
      dimension a(np, np), indx(n), vv(nmax)

      ierr = 0
      d = 1.d0
      do 12 i = 1, n
         aamax = 0.d0
         do 11 j = 1, n
            if (dabs(a(i,j)) .gt. aamax) aamax = dabs(a(i,j))
   11    continue
         if (aamax .eq. 0.d0) then
            write(unit=*,fmt=*) 'Singular matrix.'
            ierr = 1
            return
         end if
         vv(i) = 1.d0 / aamax
   12 continue

      do 19 j = 1, n
         if (j .gt. 1) then
            do 14 i = 1, j - 1
               sum = a(i,j)
               if (i .gt. 1) then
                  do 13 k = 1, i - 1
                     sum = sum - (a(i,k) * a(k,j))
   13             continue
                  a(i,j) = sum
               end if
   14       continue
         end if
         aamax = 0.d0
         do 16 i = j, n
            sum = a(i,j)
            if (j .gt. 1) then
               do 15 k = 1, j - 1
                  sum = sum - (a(i,k) * a(k,j))
   15          continue
               a(i,j) = sum
            end if
            dum = vv(i) * dabs(sum)
            if (dum .ge. aamax) then
               imax = i
               aamax = dum
            end if
   16    continue
         if (j .ne. imax) then
            do 17 k = 1, n
               dum = a(imax,k)
               a(imax,k) = a(j,k)
               a(j,k) = dum
   17       continue
            d = - d
            vv(imax) = vv(j)
         end if
         indx(j) = imax
         if (j .ne. n) then
            if (a(j,j) .eq. 0.d0) a(j,j) = tiny
               dum = 1.d0 / a(j,j)
               do 18 i = j + 1, n
                  a(i,j) = a(i,j) * dum
   18          continue
         end if
   19 continue
      if (a(n,n) .eq. 0.d0) a(n,n) = tiny
      return 
      end


      subroutine lubksb(a, n, np, indx, b)
      implicit double precision (o-z, a-h)
      dimension a(np, np), indx(n), b(n)
      ii = 0
      do 12 i = 1, n
         ll = indx(i)
         sum = b(ll)
         b(ll) = b(i)
         if (ii .ne. 0) then
            do 11 j = ii, i - 1
               sum = sum - (a(i,j) * b(j))
   11       continue
         else if (sum .ne. 0.d0) then
            ii = i
         end if
         b(i) = sum
   12 continue

      do 14 i = n, 1, -1
         sum = b(i)
         if (i .lt. n) then
            do 13 j = i + 1, n
               sum = sum - (a(i,j) * b(j))
   13       continue
         end if
         b(i) = sum / a(i,i)
   14 continue

      return 
      end