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