subroutine sgbsl(abd,lda,n,ml,mu,ipvt,b,job)
integer lda,n,ml,mu,ipvt(1),job
real abd(lda,1),b(1)
c
c sgbsl solves the real band system
c a * x = b or trans(a) * x = b
c using the factors computed by sgbco or sgbfa.
c
c on entry
c
c abd real(lda, n)
c the output from sgbco or sgbfa.
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 sgbco or sgbfa.
c
c b real(n)
c the right hand side vector.
c
c job integer
c = 0 to solve a*x = b ,
c = nonzero to solve trans(a)*x = b , where
c trans(a) is the 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 sgbco has set rcond .gt. 0.0
c or sgbfa has set info .eq. 0 .
c
c to compute inverse(a) * c where c is a matrix
c with p columns
c call sgbco(abd,lda,n,ml,mu,ipvt,rcond,z)
c if (rcond is too small) go to ...
c do 10 j = 1, p
c call sgbsl(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 saxpy,sdot
c fortran min0
c
c internal variables
c
real sdot,t
integer k,kb,l,la,lb,lm,m,nm1
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 saxpy(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 saxpy(lm,t,abd(la,k),1,b(lb),1)
40 continue
go to 100
50 continue
c
c job = nonzero, solve trans(a) * x = b
c first solve trans(u)*y = b
c
do 60 k = 1, n
lm = min0(k,m) - 1
la = m - lm
lb = k - lm
t = sdot(lm,abd(la,k),1,b(lb),1)
b(k) = (b(k) - t)/abd(m,k)
60 continue
c
c now solve trans(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) + sdot(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