subroutine csidi(a,lda,n,kpvt,det,work,job)
      integer lda,n,job
      complex a(lda,1),det(2),work(1)
      integer kpvt(1)
c
c     csidi computes the determinant and inverse
c     of a complex symmetric matrix using the factors from csifa.
c
c     on entry
c
c        a       complex(lda,n)
c                the output from csifa.
c
c        lda     integer
c                the leading dimension of the array a.
c
c        n       integer
c                the order of the matrix a.
c
c        kpvt    integer(n)
c                the pivot vector from csifa.
c
c        work    complex(n)
c                work vector.  contents destroyed.
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        a      contains the upper triangle of the inverse of
c               the original matrix.  the strict lower triangle
c               is never referenced.
c
c        det    complex(2)
c               determinant of original matrix.
c               determinant = det(1) * 10.0**det(2)
c               with 1.0 .le. abs(det(1)) .lt. 10.0
c               or det(1) = 0.0.
c
c     error condition
c
c        a division by zero may occur if the inverse is requested
c        and  csico  has set rcond .eq. 0.0
c        or  csifa  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 caxpy,ccopy,cdotu,cswap
c     fortran abs,cmplx,iabs,mod,real
c
c     internal variables.
c
      complex ak,akp1,akkp1,cdotu,d,t,temp
      real ten
      integer j,jb,k,km1,ks,kstep
      logical noinv,nodet
c
      complex zdum
      real cabs1
      cabs1(zdum) = abs(real(zdum)) + abs(aimag(zdum))
c
      noinv = mod(job,10) .eq. 0
      nodet = mod(job,100)/10 .eq. 0
c
      if (nodet) go to 100
         det(1) = (1.0e0,0.0e0)
         det(2) = (0.0e0,0.0e0)
         ten = 10.0e0
         t = (0.0e0,0.0e0)
         do 90 k = 1, n
            d = a(k,k)
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.0e0) go to 10
                  t = a(k,k+1)
                  d = (d/t)*a(k+1,k+1) - t
               go to 20
   10          continue
                  d = t
                  t = (0.0e0,0.0e0)
   20          continue
   30       continue
c
            det(1) = d*det(1)
            if (cabs1(det(1)) .eq. 0.0e0) go to 80
   40          if (cabs1(det(1)) .ge. 1.0e0) go to 50
                  det(1) = cmplx(ten,0.0e0)*det(1)
                  det(2) = det(2) - (1.0e0,0.0e0)
               go to 40
   50          continue
   60          if (cabs1(det(1)) .lt. ten) go to 70
                  det(1) = det(1)/cmplx(ten,0.0e0)
                  det(2) = det(2) + (1.0e0,0.0e0)
               go to 60
   70          continue
   80       continue
   90    continue
  100 continue
c
c     compute inverse(a)
c
      if (noinv) go to 230
         k = 1
  110    if (k .gt. n) go to 220
            km1 = k - 1
            if (kpvt(k) .lt. 0) go to 140
c
c              1 by 1
c
               a(k,k) = (1.0e0,0.0e0)/a(k,k)
               if (km1 .lt. 1) go to 130
                  call ccopy(km1,a(1,k),1,work,1)
                  do 120 j = 1, km1
                     a(j,k) = cdotu(j,a(1,j),1,work,1)
                     call caxpy(j-1,work(j),a(1,j),1,a(1,k),1)
  120             continue
                  a(k,k) = a(k,k) + cdotu(km1,work,1,a(1,k),1)
  130          continue
               kstep = 1
            go to 180
  140       continue
c
c              2 by 2
c
               t = a(k,k+1)
               ak = a(k,k)/t
               akp1 = a(k+1,k+1)/t
               akkp1 = a(k,k+1)/t
               d = t*(ak*akp1 - (1.0e0,0.0e0))
               a(k,k) = akp1/d
               a(k+1,k+1) = ak/d
               a(k,k+1) = -akkp1/d
               if (km1 .lt. 1) go to 170
                  call ccopy(km1,a(1,k+1),1,work,1)
                  do 150 j = 1, km1
                     a(j,k+1) = cdotu(j,a(1,j),1,work,1)
                     call caxpy(j-1,work(j),a(1,j),1,a(1,k+1),1)
  150             continue
                  a(k+1,k+1) = a(k+1,k+1)
     *                         + cdotu(km1,work,1,a(1,k+1),1)
                  a(k,k+1) = a(k,k+1) + cdotu(km1,a(1,k),1,a(1,k+1),1)
                  call ccopy(km1,a(1,k),1,work,1)
                  do 160 j = 1, km1
                     a(j,k) = cdotu(j,a(1,j),1,work,1)
                     call caxpy(j-1,work(j),a(1,j),1,a(1,k),1)
  160             continue
                  a(k,k) = a(k,k) + cdotu(km1,work,1,a(1,k),1)
  170          continue
               kstep = 2
  180       continue
c
c           swap
c
            ks = iabs(kpvt(k))
            if (ks .eq. k) go to 210
               call cswap(ks,a(1,ks),1,a(1,k),1)
               do 190 jb = ks, k
                  j = k + ks - jb
                  temp = a(j,k)
                  a(j,k) = a(ks,j)
                  a(ks,j) = temp
  190          continue
               if (kstep .eq. 1) go to 200
                  temp = a(ks,k+1)
                  a(ks,k+1) = a(k,k+1)
                  a(k,k+1) = temp
  200          continue
  210       continue
            k = k + kstep
         go to 110
  220    continue
  230 continue
      return
      end