subroutine zgbsl(abd,lda,n,ml,mu,ipvt,b,job) integer lda,n,ml,mu,ipvt(1),job complex*16 abd(lda,1),b(1) c c zgbsl solves the complex*16 band system c a * x = b or ctrans(a) * x = b c using the factors computed by zgbco or zgbfa. c c on entry c c abd complex*16(lda, n) c the output from zgbco or zgbfa. c c lda integer c the leading dimension of the array abd . c c n integer c the order of the original matrix. c c ml integer c number of diagonals below the main diagonal. c c mu integer c number of diagonals above the main diagonal. c c ipvt integer(n) c the pivot vector from zgbco or zgbfa. 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 zgbco has set rcond .gt. 0.0 c or zgbfa has set info .eq. 0 . c c to compute inverse(a) * c where c is a matrix c with p columns c call zgbco(abd,lda,n,ml,mu,ipvt,rcond,z) c if (rcond is too small) go to ... c do 10 j = 1, p c call zgbsl(abd,lda,n,ml,mu,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,min0 c c internal variables c complex*16 zdotc,t integer k,kb,l,la,lb,lm,m,nm1 double precision dreal,dimag complex*16 zdumr,zdumi dreal(zdumr) = zdumr dimag(zdumi) = (0.0d0,-1.0d0)*zdumi c m = mu + ml + 1 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 (ml .eq. 0) go to 30 if (nm1 .lt. 1) go to 30 do 20 k = 1, nm1 lm = min0(ml,n-k) l = ipvt(k) t = b(l) if (l .eq. k) go to 10 b(l) = b(k) b(k) = t 10 continue call zaxpy(lm,t,abd(m+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)/abd(m,k) lm = min0(k,m) - 1 la = m - lm lb = k - lm t = -b(k) call zaxpy(lm,t,abd(la,k),1,b(lb),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 lm = min0(k,m) - 1 la = m - lm lb = k - lm t = zdotc(lm,abd(la,k),1,b(lb),1) b(k) = (b(k) - t)/dconjg(abd(m,k)) 60 continue c c now solve ctrans(l)*x = y c if (ml .eq. 0) go to 90 if (nm1 .lt. 1) go to 90 do 80 kb = 1, nm1 k = n - kb lm = min0(ml,n-k) b(k) = b(k) + zdotc(lm,abd(m+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