subroutine LUDCMP(a,n,np,indx,d)
      implicit double precision (a-h,o-z)
c
c Given an NxN matrix A, with dimension NP, this routine replaces it by
c the LU decomposition ofa rowwise permutation of itself.  A and N are
c input.  A is the output, arranged in equation 2.3.14 on p. 34.  INDX is
c an output vector which records the row permutation resulting from
c partial pivoting.  D is output as +/- 1 depending on whether the number
c of row interchanges was even or odd, respectively.  Used in combination
c with LUBKSB to solve linear equations or invert a matrix.
c 
c Numerical Recipes, p. 38-9.
c
      INTEGER n,np,indx(n),NMAX,i,imax,j,k
      DOUBLE PRECISION d,a(np,np),TINY
      PARAMETER(nmax=500, tiny=1.d-40)
      DOUBLE PRECISION aamax,dum,sum,vv(NMAX)
      d=1.0
      do 12 i=1,n
        aamax=0.d0
        do 11 j=1,n
          if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j)) 
11      continue
        if (aamax.eq.0.d0) PAUSE ' Singular matrix in LUDCMP'
        vv(i)=1.d0/aamax
12    continue
      do 19 j=1,n
        do 14 i=1,j-1
          sum=a(i,j)
          do 13 k=1,i-1
            sum=sum-a(i,k)*a(k,j)
13        continue
          a(i,j)=sum
14      continue
        aamax=0.d0
        do 16 i=j,n
          sum=a(i,j)
          do 15 k=1,j-1
            sum=sum-a(i,k)*a(k,j)
15        continue
          a(i,j)=sum
          dum=vv(i)*abs(sum)
          if (dum.ge.aamax) then
            imax=i
            aamax=dum
          endif
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)
        endif
        indx(j)=imax
        if (a(j,j).eq.0.d0) a(j,j)=tiny
        if (j.ne.n) then
          dum=1.d0/a(j,j)
          do 18 i=j+1,n
            a(i,j)=a(i,j)*dum
18        continue
        endif
19    continue
      return
      end
c
      subroutine LUBKSB(a,n,np,indx,b)
      implicit double precision (a-h,o-z)
c
c Solves the set of N linear equations A.X=B. Here A is input, not as
c the matrix A but rather its LU decomposition, determined by the routine
c LUDCMP.  INDX is input as the permutation vector returned by LUDCMP.  B
c is input as teh right hand side vector B and returns with the solution
c vector X.  A, N, NP, and INDX are not modified by this routine and can
c be left in place for successive calls with different right-hand sides B.
c This routie takes into account the possibility that B will begin with
c many zero elements, so it is efficient for use in matrix inversion.
c
c Numerical Recipes, p. 39
c
      INTEGER n,np,indx(n)
      DOUBLE PRECISION a(np,np),b(n)
      INTEGER i,ii,j,ll
      DOUBLE PRECISION sum
      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
        endif
        b(i)=sum
12    continue
      do 14 i=n,1,-1
        sum=b(i)
        do 13 j=i+1,n
          sum=sum-a(i,j)*b(j)
13      continue
        b(i)=sum/a(i,i)
14    continue
      return
      end