subroutine dpbfa(abd,lda,n,m,info)
integer lda,n,m,info
double precision abd(lda,1)
c
c dpbfa factors a double precision symmetric positive definite
c matrix stored in band form.
c
c dpbfa is usually called by dpbco, 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 the matrix to be factored. the columns of the upper
c triangle are stored in the columns of abd and the
c diagonals of the upper triangle are stored in the
c rows of abd . see the comments below for details.
c
c lda integer
c the leading dimension of the array abd .
c lda must be .ge. m + 1 .
c
c n integer
c the order of the matrix a .
c
c m integer
c the number of diagonals above the main diagonal.
c 0 .le. m .lt. n .
c
c on return
c
c abd an upper triangular matrix r , stored in band
c form, so that a = trans(r)*r .
c
c info integer
c = 0 for normal return.
c = k if the leading minor of order k is not
c positive definite.
c
c band storage
c
c if a is a symmetric positive definite band matrix,
c the following program segment will set up the input.
c
c m = (band width above diagonal)
c do 20 j = 1, n
c i1 = max0(1, j-m)
c do 10 i = i1, j
c k = i-j+m+1
c abd(k,j) = a(i,j)
c 10 continue
c 20 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 ddot
c fortran max0,dsqrt
c
c internal variables
c
double precision ddot,t
double precision s
integer ik,j,jk,k,mu
c begin block with ...exits to 40
c
c
do 30 j = 1, n
info = j
s = 0.0d0
ik = m + 1
jk = max0(j-m,1)
mu = max0(m+2-j,1)
if (m .lt. mu) go to 20
do 10 k = mu, m
t = abd(k,j) - ddot(k-mu,abd(ik,jk),1,abd(mu,j),1)
t = t/abd(m+1,jk)
abd(k,j) = t
s = s + t*t
ik = ik - 1
jk = jk + 1
10 continue
20 continue
s = abd(m+1,j) - s
c ......exit
if (s .le. 0.0d0) go to 40
abd(m+1,j) = dsqrt(s)
30 continue
info = 0
40 continue
return
end