subroutine zgesl(a,lda,n,ipvt,b,job) integer lda,n,ipvt(1),job complex*16 a(lda,1),b(1) c c zgesl solves the complex*16 system c a * x = b or ctrans(a) * x = b 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 b complex*16(n) c the right hand side vector. c c job integer c = 0 to solve a*x = b , c = nonzero to solve ctrans(a)*x = b where c ctrans(a) is the conjugate transpose. c c on return c c b the solution vector x . c c error condition c c a division by zero will occur if the input factor contains a c zero on the diagonal. technically this indicates singularity c but it is often caused by improper arguments or improper c setting of lda . it will not occur if the subroutines are c called correctly and if zgeco has set rcond .gt. 0.0 c or zgefa has set info .eq. 0 . c c to compute inverse(a) * c where c is a matrix c with p columns c call zgeco(a,lda,n,ipvt,rcond,z) c if (rcond is too small) go to ... c do 10 j = 1, p c call zgesl(a,lda,n,ipvt,c(1,j),0) c 10 continue 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,zdotc c fortran dconjg c c internal variables c complex*16 zdotc,t integer k,kb,l,nm1 double precision dreal,dimag complex*16 zdumr,zdumi dreal(zdumr) = zdumr dimag(zdumi) = (0.0d0,-1.0d0)*zdumi c nm1 = n - 1 if (job .ne. 0) go to 50 c c job = 0 , solve a * x = b c first solve l*y = b c if (nm1 .lt. 1) go to 30 do 20 k = 1, nm1 l = ipvt(k) t = b(l) if (l .eq. k) go to 10 b(l) = b(k) b(k) = t 10 continue call zaxpy(n-k,t,a(k+1,k),1,b(k+1),1) 20 continue 30 continue c c now solve u*x = y c do 40 kb = 1, n k = n + 1 - kb b(k) = b(k)/a(k,k) t = -b(k) call zaxpy(k-1,t,a(1,k),1,b(1),1) 40 continue go to 100 50 continue c c job = nonzero, solve ctrans(a) * x = b c first solve ctrans(u)*y = b c do 60 k = 1, n t = zdotc(k-1,a(1,k),1,b(1),1) b(k) = (b(k) - t)/dconjg(a(k,k)) 60 continue c c now solve ctrans(l)*x = y c if (nm1 .lt. 1) go to 90 do 80 kb = 1, nm1 k = n - kb b(k) = b(k) + zdotc(n-k,a(k+1,k),1,b(k+1),1) l = ipvt(k) if (l .eq. k) go to 70 t = b(l) b(l) = b(k) b(k) = t 70 continue 80 continue 90 continue 100 continue return end