subroutine zspdi(ap,n,kpvt,det,work,job)
      integer n,job
      complex*16 ap(1),work(1),det(2)
      integer kpvt(1)
c
c     zspdi computes the determinant and inverse
c     of a complex*16 symmetric matrix using the factors from zspfa,
c     where the matrix is stored in packed form.
c
c     on entry
c
c        ap      complex*16 (n*(n+1)/2)
c                the output from zspfa.
c
c        n       integer
c                the order of the matrix a.
c
c        kpvt    integer(n)
c                the pivot vector from zspfa.
c
c        work    complex*16(n)
c                work vector.  contents ignored.
c
c        job     integer
c                job has the decimal expansion  ab  where
c                   if  b .ne. 0, the inverse is computed,
c                   if  a .ne. 0, the determinant is computed,
c
c                for example, job = 11  gives both.
c
c     on return
c
c        variables not requested by job are not used.
c
c        ap     contains the upper triangle of the inverse of
c               the original matrix, stored in packed form.
c               the columns of the upper triangle are stored
c               sequentially in a one-dimensional array.
c
c        det    complex*16(2)
c               determinant of original matrix.
c               determinant = det(1) * 10.0**det(2)
c               with 1.0 .le. dabs(det(1)) .lt. 10.0
c               or det(1) = 0.0.
c
c     error condition
c
c        a division by zero will occur if the inverse is requested
c        and  zspco  has set rcond .eq. 0.0
c        or  zspfa  has set  info .ne. 0 .
c
c     linpack. this version dated 08/14/78 .
c     james bunch, univ. calif. san diego, argonne nat. lab.
c
c     subroutines and functions
c
c     blas zaxpy,zcopy,zdotu,zswap
c     fortran dabs,dcmplx,iabs,mod
c
c     internal variables.
c
      complex*16 ak,akkp1,akp1,zdotu,d,t,temp
      double precision ten
      integer ij,ik,ikp1,iks,j,jb,jk,jkp1
      integer k,kk,kkp1,km1,ks,ksj,kskp1,kstep
      logical noinv,nodet
c
      complex*16 zdum
      double precision cabs1
      double precision dreal,dimag
      complex*16 zdumr,zdumi
      dreal(zdumr) = zdumr
      dimag(zdumi) = (0.0d0,-1.0d0)*zdumi
      cabs1(zdum) = dabs(dreal(zdum)) + dabs(dimag(zdum))
c
      noinv = mod(job,10) .eq. 0
      nodet = mod(job,100)/10 .eq. 0
c
      if (nodet) go to 110
         det(1) = (1.0d0,0.0d0)
         det(2) = (0.0d0,0.0d0)
         ten = 10.0d0
         t = (0.0d0,0.0d0)
         ik = 0
         do 100 k = 1, n
            kk = ik + k
            d = ap(kk)
c
c           check if 1 by 1
c
            if (kpvt(k) .gt. 0) go to 30
c
c              2 by 2 block
c              use det (d  t)  =  (d/t * c - t) * t
c                      (t  c)
c              to avoid underflow/overflow troubles.
c              take two passes through scaling.  use  t  for flag.
c
               if (cabs1(t) .ne. 0.0d0) go to 10
                  ikp1 = ik + k
                  kkp1 = ikp1 + k
                  t = ap(kkp1)
                  d = (d/t)*ap(kkp1+1) - t
               go to 20
   10          continue
                  d = t
                  t = (0.0d0,0.0d0)
   20          continue
   30       continue
c
            if (nodet) go to 90
               det(1) = d*det(1)
               if (cabs1(det(1)) .eq. 0.0d0) go to 80
   40             if (cabs1(det(1)) .ge. 1.0d0) go to 50
                     det(1) = dcmplx(ten,0.0d0)*det(1)
                     det(2) = det(2) - (1.0d0,0.0d0)
                  go to 40
   50             continue
   60             if (cabs1(det(1)) .lt. ten) go to 70
                     det(1) = det(1)/dcmplx(ten,0.0d0)
                     det(2) = det(2) + (1.0d0,0.0d0)
                  go to 60
   70             continue
   80          continue
   90       continue
            ik = ik + k
  100    continue
  110 continue
c
c     compute inverse(a)
c
      if (noinv) go to 240
         k = 1
         ik = 0
  120    if (k .gt. n) go to 230
            km1 = k - 1
            kk = ik + k
            ikp1 = ik + k
            if (kpvt(k) .lt. 0) go to 150
c
c              1 by 1
c
               ap(kk) = (1.0d0,0.0d0)/ap(kk)
               if (km1 .lt. 1) go to 140
                  call zcopy(km1,ap(ik+1),1,work,1)
                  ij = 0
                  do 130 j = 1, km1
                     jk = ik + j
                     ap(jk) = zdotu(j,ap(ij+1),1,work,1)
                     call zaxpy(j-1,work(j),ap(ij+1),1,ap(ik+1),1)
                     ij = ij + j
  130             continue
                  ap(kk) = ap(kk) + zdotu(km1,work,1,ap(ik+1),1)
  140          continue
               kstep = 1
            go to 190
  150       continue
c
c              2 by 2
c
               kkp1 = ikp1 + k
               t = ap(kkp1)
               ak = ap(kk)/t
               akp1 = ap(kkp1+1)/t
               akkp1 = ap(kkp1)/t
               d = t*(ak*akp1 - (1.0d0,0.0d0))
               ap(kk) = akp1/d
               ap(kkp1+1) = ak/d
               ap(kkp1) = -akkp1/d
               if (km1 .lt. 1) go to 180
                  call zcopy(km1,ap(ikp1+1),1,work,1)
                  ij = 0
                  do 160 j = 1, km1
                     jkp1 = ikp1 + j
                     ap(jkp1) = zdotu(j,ap(ij+1),1,work,1)
                     call zaxpy(j-1,work(j),ap(ij+1),1,ap(ikp1+1),1)
                     ij = ij + j
  160             continue
                  ap(kkp1+1) = ap(kkp1+1)
     *                         + zdotu(km1,work,1,ap(ikp1+1),1)
                  ap(kkp1) = ap(kkp1)
     *                       + zdotu(km1,ap(ik+1),1,ap(ikp1+1),1)
                  call zcopy(km1,ap(ik+1),1,work,1)
                  ij = 0
                  do 170 j = 1, km1
                     jk = ik + j
                     ap(jk) = zdotu(j,ap(ij+1),1,work,1)
                     call zaxpy(j-1,work(j),ap(ij+1),1,ap(ik+1),1)
                     ij = ij + j
  170             continue
                  ap(kk) = ap(kk) + zdotu(km1,work,1,ap(ik+1),1)
  180          continue
               kstep = 2
  190       continue
c
c           swap
c
            ks = iabs(kpvt(k))
            if (ks .eq. k) go to 220
               iks = (ks*(ks - 1))/2
               call zswap(ks,ap(iks+1),1,ap(ik+1),1)
               ksj = ik + ks
               do 200 jb = ks, k
                  j = k + ks - jb
                  jk = ik + j
                  temp = ap(jk)
                  ap(jk) = ap(ksj)
                  ap(ksj) = temp
                  ksj = ksj - (j - 1)
  200          continue
               if (kstep .eq. 1) go to 210
                  kskp1 = ikp1 + ks
                  temp = ap(kskp1)
                  ap(kskp1) = ap(kkp1)
                  ap(kkp1) = temp
  210          continue
  220       continue
            ik = ik + k
            if (kstep .eq. 2) ik = ik + k + 1
            k = k + kstep
         go to 120
  230    continue
  240 continue
      return
      end