dpbfa.f

      subroutine dpbfa(abd,lda,n,m,info)

      implicit none
      integer lda,n,m,info
      double precision abd(lda,n)
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,sqrt
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.0e0
	    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.0e0) go to 40

	    abd(m+1,j) = sqrt(s)

   30	 continue
	 info = 0
   40 continue
      return
      end