subroutine ludcmp(a,n,np,indx,d)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      integer n,np,indx(n),NMAX
      double precision d, a(np,np), TINY
      parameter (NMAX=500,TINY=1d-20)
      integer i,imax,j,k
      double precision aamax,dum,sum,vv(NMAX)

      
      d=1.
      do i=1,n
         aamax=0.
         do j=1,n
            if(dabs(a(i,j)).gt.aamax) aamax=abs(a(i,j))
         end do
         if(aamax.eq.0.0d0) then
            print *,'singular matrix in ludcmp'
            return
         endif
         vv(i) = 1.d0/aamax
      end do

      do j=1,n
         do i=1,j-1
            sum=a(i,j)
            do k=1,i-1
               sum=sum-a(i,k)*a(k,j)
            end do
            a(i,j)=sum
         end do
         aamax=0.
         do i=j,n
            sum=a(i,j)
            do k=1,j-1
               sum=sum-a(i,k)*a(k,j)
            end do
            a(i,j)=sum
            dum=vv(i)*abs(sum)
            if(dum.ge.aamax) then
               imax=i
               aamax=dum
            end if
         end do
         if(j.ne.imax) then
            do k=1,n
               dum=a(imax,k)
               a(imax,k)=a(j,k)
               a(j,k)=dum
            end do
            d=-d
            vv(imax)=vv(j)
         end if
         indx(j)=imax
         if(a(j,j).eq.0.0d0) a(j,j)=TINY
         if(j.ne.n)then
            dum=1./a(j,j)
            do i=j+1,n
               a(i,j)=a(i,j)*dum
            end do
         end if
      end do
      return
      end


      subroutine lubksb(a,n,np,indx,b)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      integer n,np,indx(n)
      double precision a(np,np),b(n)

      integer i,ii,j,ll
      double precision sum
      
      ii=0
      do i=1,n
         ll=indx(i)
         sum=b(ll)
         b(ll)=b(i)
         if(ii.ne.0) then
            do j=ii,i-1
               sum=sum-a(i,j)*b(j)
            end do
         else if (sum.ne.0.0d0) then
            ii=i
         end if
         b(i)=sum
      end do
      do i=n,1,-1
         sum=b(i)
         do j=i+1,n
            sum=sum-a(i,j)*b(j)
         end do
         b(i)=sum/a(i,i)
      end do
      return
      end