subroutine zgedi(a,lda,n,ipvt,det,work,job) integer lda,n,ipvt(1),job complex*16 a(lda,1),det(2),work(1) c c zgedi computes the determinant and inverse of a matrix c using the factors computed by zgeco or zgefa. c c on entry c c a complex*16(lda, n) c the output from zgeco or zgefa. 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 ipvt integer(n) c the pivot vector from zgeco or zgefa. c c work complex*16(n) c work vector. contents destroyed. c c job integer c = 11 both determinant and inverse. c = 01 inverse only. c = 10 determinant only. c c on return c c a inverse of original matrix if requested. c otherwise unchanged. c c det complex*16(2) c determinant of original matrix if requested. c otherwise not referenced. c determinant = det(1) * 10.0**det(2) c with 1.0 .le. cabs1(det(1)) .lt. 10.0 c or det(1) .eq. 0.0 . c c error condition c c a division by zero will occur if the input factor contains c a zero on the diagonal and the inverse is requested. c it will not occur if the subroutines are called correctly c and if zgeco has set rcond .gt. 0.0 or zgefa has set c info .eq. 0 . c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas zaxpy,zscal,zswap c fortran dabs,dcmplx,mod c c internal variables c complex*16 t double precision ten integer i,j,k,kb,kp1,l,nm1 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 c compute determinant c if (job/10 .eq. 0) go to 70 det(1) = (1.0d0,0.0d0) det(2) = (0.0d0,0.0d0) ten = 10.0d0 do 50 i = 1, n if (ipvt(i) .ne. i) det(1) = -det(1) det(1) = a(i,i)*det(1) c ...exit if (cabs1(det(1)) .eq. 0.0d0) go to 60 10 if (cabs1(det(1)) .ge. 1.0d0) go to 20 det(1) = dcmplx(ten,0.0d0)*det(1) det(2) = det(2) - (1.0d0,0.0d0) go to 10 20 continue 30 if (cabs1(det(1)) .lt. ten) go to 40 det(1) = det(1)/dcmplx(ten,0.0d0) det(2) = det(2) + (1.0d0,0.0d0) go to 30 40 continue 50 continue 60 continue 70 continue c c compute inverse(u) c if (mod(job,10) .eq. 0) go to 150 do 100 k = 1, n a(k,k) = (1.0d0,0.0d0)/a(k,k) t = -a(k,k) call zscal(k-1,t,a(1,k),1) kp1 = k + 1 if (n .lt. kp1) go to 90 do 80 j = kp1, n t = a(k,j) a(k,j) = (0.0d0,0.0d0) call zaxpy(k,t,a(1,k),1,a(1,j),1) 80 continue 90 continue 100 continue c c form inverse(u)*inverse(l) c nm1 = n - 1 if (nm1 .lt. 1) go to 140 do 130 kb = 1, nm1 k = n - kb kp1 = k + 1 do 110 i = kp1, n work(i) = a(i,k) a(i,k) = (0.0d0,0.0d0) 110 continue do 120 j = kp1, n t = work(j) call zaxpy(n,t,a(1,j),1,a(1,k),1) 120 continue l = ipvt(k) if (l .ne. k) call zswap(n,a(1,k),1,a(1,l),1) 130 continue 140 continue 150 continue return end