subroutine dgbfa(abd,lda,n,ml,mu,ipvt,info)
integer lda,n,ml,mu,ipvt(1),info
double precision abd(lda,1)
c
c dgbfa factors a double precision band matrix by elimination.
c
c dgbfa is usually called by dgbco, but it can be called
c directly with a saving in time if rcond is not needed.
c
c on entry
c
c abd double precision(lda, n)
c contains the matrix in band storage. the columns
c of the matrix are stored in the columns of abd and
c the diagonals of the matrix are stored in rows
c ml+1 through 2*ml+mu+1 of abd .
c see the comments below for details.
c
c lda integer
c the leading dimension of the array abd .
c lda must be .ge. 2*ml + mu + 1 .
c
c n integer
c the order of the original matrix.
c
c ml integer
c number of diagonals below the main diagonal.
c 0 .le. ml .lt. n .
c
c mu integer
c number of diagonals above the main diagonal.
c 0 .le. mu .lt. n .
c more efficient if ml .le. mu .
c on return
c
c abd an upper triangular matrix in band storage and
c the multipliers which were used to obtain it.
c the factorization can be written a = l*u where
c l is a product of permutation and unit lower
c triangular matrices and u is upper triangular.
c
c ipvt integer(n)
c an integer vector of pivot indices.
c
c info integer
c = 0 normal value.
c = k if u(k,k) .eq. 0.0 . this is not an error
c condition for this subroutine, but it does
c indicate that dgbsl will divide by zero if
c called. use rcond in dgbco for a reliable
c indication of singularity.
c
c band storage
c
c if a is a band matrix, the following program segment
c will set up the input.
c
c ml = (band width below the diagonal)
c mu = (band width above the diagonal)
c m = ml + mu + 1
c do 20 j = 1, n
c i1 = max0(1, j-mu)
c i2 = min0(n, j+ml)
c do 10 i = i1, i2
c k = i - j + m
c abd(k,j) = a(i,j)
c 10 continue
c 20 continue
c
c this uses rows ml+1 through 2*ml+mu+1 of abd .
c in addition, the first ml rows in abd are used for
c elements generated during the triangularization.
c the total number of rows needed in abd is 2*ml+mu+1 .
c the ml+mu by ml+mu upper left triangle and the
c ml by ml lower right triangle are not referenced.
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 daxpy,dscal,idamax
c fortran max0,min0
c
c internal variables
c
double precision t
integer i,idamax,i0,j,ju,jz,j0,j1,k,kp1,l,lm,m,mm,nm1
c
c
m = ml + mu + 1
info = 0
c
c zero initial fill-in columns
c
j0 = mu + 2
j1 = min0(n,m) - 1
if (j1 .lt. j0) go to 30
do 20 jz = j0, j1
i0 = m + 1 - jz
do 10 i = i0, ml
abd(i,jz) = 0.0d0
10 continue
20 continue
30 continue
jz = j1
ju = 0
c
c gaussian elimination with partial pivoting
c
nm1 = n - 1
if (nm1 .lt. 1) go to 130
do 120 k = 1, nm1
kp1 = k + 1
c
c zero next fill-in column
c
jz = jz + 1
if (jz .gt. n) go to 50
if (ml .lt. 1) go to 50
do 40 i = 1, ml
abd(i,jz) = 0.0d0
40 continue
50 continue
c
c find l = pivot index
c
lm = min0(ml,n-k)
l = idamax(lm+1,abd(m,k),1) + m - 1
ipvt(k) = l + k - m
c
c zero pivot implies this column already triangularized
c
if (abd(l,k) .eq. 0.0d0) go to 100
c
c interchange if necessary
c
if (l .eq. m) go to 60
t = abd(l,k)
abd(l,k) = abd(m,k)
abd(m,k) = t
60 continue
c
c compute multipliers
c
t = -1.0d0/abd(m,k)
call dscal(lm,t,abd(m+1,k),1)
c
c row elimination with column indexing
c
ju = min0(max0(ju,mu+ipvt(k)),n)
mm = m
if (ju .lt. kp1) go to 90
do 80 j = kp1, ju
l = l - 1
mm = mm - 1
t = abd(l,j)
if (l .eq. mm) go to 70
abd(l,j) = abd(mm,j)
abd(mm,j) = t
70 continue
call daxpy(lm,t,abd(m+1,k),1,abd(mm+1,j),1)
80 continue
90 continue
go to 110
100 continue
info = k
110 continue
120 continue
130 continue
ipvt(n) = n
if (abd(m,n) .eq. 0.0d0) info = n
return
end