subroutine cgesl(a,lda,n,ipvt,b,job)
integer lda,n,ipvt(1),job
complex a(lda,1),b(1)
c
c cgesl solves the complex system
c a * x = b or ctrans(a) * x = b
c using the factors computed by cgeco or cgefa.
c
c on entry
c
c a complex(lda, n)
c the output from cgeco or cgefa.
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 cgeco or cgefa.
c
c b complex(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 cgeco has set rcond .gt. 0.0
c or cgefa has set info .eq. 0 .
c
c to compute inverse(a) * c where c is a matrix
c with p columns
c call cgeco(a,lda,n,ipvt,rcond,z)
c if (rcond is too small) go to ...
c do 10 j = 1, p
c call cgesl(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 caxpy,cdotc
c fortran conjg
c
c internal variables
c
complex cdotc,t
integer k,kb,l,nm1
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 caxpy(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 caxpy(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 = cdotc(k-1,a(1,k),1,b(1),1)
b(k) = (b(k) - t)/conjg(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) + cdotc(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