subroutine ztrdi(t,ldt,n,det,job,info) integer ldt,n,job,info complex*16 t(ldt,1),det(2) c c ztrdi computes the determinant and inverse of a complex*16 c triangular matrix. c c on entry c c t complex*16(ldt,n) c t contains the triangular matrix. the zero c elements of the matrix are not referenced, and c the corresponding elements of the array can be c used to store other information. c c ldt integer c ldt is the leading dimension of the array t. c c n integer c n is the order of the system. c c job integer c = 010 no det, inverse of lower triangular. c = 011 no det, inverse of upper triangular. c = 100 det, no inverse. c = 110 det, inverse of lower triangular. c = 111 det, inverse of upper triangular. c c on return c c t 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 info integer c info contains zero if the system is nonsingular c and the inverse is requested. c otherwise info contains the index of c a zero diagonal element of t. c 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 c fortran dabs,dcmplx,mod c c internal variables c complex*16 temp double precision ten integer i,j,k,kb,km1,kp1 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 begin block permitting ...exits to 180 c c compute determinant c if (job/100 .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 det(1) = t(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 of upper triangular c if (mod(job/10,10) .eq. 0) go to 170 if (mod(job,10) .eq. 0) go to 120 c begin block permitting ...exits to 110 do 100 k = 1, n info = k c ......exit if (cabs1(t(k,k)) .eq. 0.0d0) go to 110 t(k,k) = (1.0d0,0.0d0)/t(k,k) temp = -t(k,k) call zscal(k-1,temp,t(1,k),1) kp1 = k + 1 if (n .lt. kp1) go to 90 do 80 j = kp1, n temp = t(k,j) t(k,j) = (0.0d0,0.0d0) call zaxpy(k,temp,t(1,k),1,t(1,j),1) 80 continue 90 continue 100 continue info = 0 110 continue go to 160 120 continue c c compute inverse of lower triangular c do 150 kb = 1, n k = n + 1 - kb info = k c ............exit if (cabs1(t(k,k)) .eq. 0.0d0) go to 180 t(k,k) = (1.0d0,0.0d0)/t(k,k) temp = -t(k,k) if (k .ne. n) call zscal(n-k,temp,t(k+1,k),1) km1 = k - 1 if (km1 .lt. 1) go to 140 do 130 j = 1, km1 temp = t(k,j) t(k,j) = (0.0d0,0.0d0) call zaxpy(n-k+1,temp,t(k,k),1,t(k,j),1) 130 continue 140 continue 150 continue info = 0 160 continue 170 continue 180 continue return end