      subroutine cludcmp(a,n,np,indx,d)
c**********************************************************************
c**** Perform LU decomposition on complex matrix A                 ****
c**** Input:  a(np,np)                                             ****
c****         np is physical dimension of A                        ****
c****         n actual dimension                                   ****
c**** Output: indx(n) is a vector recording row permutations       ****
c****         d = +1/-1  even/odd interchanges                     ****
c**********************************************************************
      integer          d,i,j,k,imax,n,np,indx(np)
      double complex   a(np,np)
      double precision vv(np)
      double complex   sum
      double precision aamax,tiny,dum
      intrinsic        cdabs

      tiny=1.0d-20
      d = 1
      do 12 i=1,n
        aamax=0.0d0
        do 11 j=1,n
          if (cdabs(a(i,j)).gt.aamax) aamax = cdabs(a(i,j))
  11    continue
        if (aamax.eq.0.0d0) pause 'Singular matrix'
        vv(i) = 1.0d0/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.0d0
         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)*cdabs(sum)
           if (dum.ge.aamax) then
             imax = i
             aamax = dum
           endif
  16     continue
         if (j.ne.imax) then
           do 17 k=1,n
             sum = a(imax,k)
             a(imax,k) = a(j,k)
             a(j,k) = sum
  17       continue
           d = -d
           vv(imax) = vv(j)
         endif
         indx(j) = imax
         if (j.ne.n) then
           if (cdabs(a(j,j)).eq.0.0d0) a(j,j) = tiny
           sum = 1.0d0/a(j,j)
           do 18 i=j+1,n
             a(i,j) = a(i,j)*sum
  18       continue
         endif
  19   continue
       if (cdabs(a(n,n)).eq.0.0d0) a(n,n) = tiny
       return
       end
