subroutine bchfac ( w, nbands, nrow, diag )
c  from  * a practical guide to splines *  by c. de boor    
constructs cholesky factorization
c                     c  =  l * d * l-transpose
c  with l unit lower triangular and d diagonal, for given matrix c of
c  order  n r o w , in case  c  is (symmetric) positive semidefinite
c  and  b a n d e d , having  n b a n d s  diagonals at and below the
c  main diagonal.
c
c******  i n p u t  ******
c  nrow.....is the order of the matrix  c .
c  nbands.....indicates its bandwidth, i.e.,
c          c(i,j) = 0 for i-j .ge. nbands .
c  w.....workarray of size (nbands,nrow)  containing the  nbands  diago-
c        nals in its rows, with the main diagonal in row  1 . precisely,
c        w(i,j)  contains  c(i+j-1,j), i=1,...,nbands, j=1,...,nrow.
c          for example, the interesting entries of a seven diagonal sym-
c        metric matrix  c  of order  9  would be stored in  w  as
c
c
c
c
c
c
c        all other entries of  w  not identified in this way with an en-
c        try of  c  are never referenced .
c  diag.....is a work array of length  nrow .
c
c******  o u t p u t  ******
c  w.....contains the cholesky factorization  c = l*d*l-transp, with
c        w(1,i) containing  1/d(i,i)
c        and  w(i,j)  containing  l(i-1+j,j), i=2,...,nbands.
c
c******  m e t h o d  ******
c   gauss elimination, adapted to the symmetry and bandedness of  c , is
c   used .
c     near zero pivots are handled in a special way. the diagonal ele-
c  ment c(n,n) = w(1,n) is saved initially in  diag(n), all n. at the n-
c  th elimination step, the current pivot element, viz.  w(1,n), is com-
c  pared with its original value, diag(n). if, as the result of prior
c  elimination steps, this element has been reduced by about a word
c  length, (i.e., if w(1,n)+diag(n) .le. diag(n)), then the pivot is de-
c  clared to be zero, and the entire n-th row is declared to be linearly
c  dependent on the preceding rows. this has the effect of producing
c   x(n) = 0  when solving  c*x = b  for  x, regardless of  b. justific-
c  ation for this is as follows. in contemplated applications of this
c  program, the given equations are the normal equations for some least-
c  squares approximation problem, diag(n) = c(n,n) gives the norm-square
c  of the n-th basis function, and, at this point,  w(1,n)  contains the
c  norm-square of the error in the least-squares approximation to the n-
c  th basis function by linear combinations of the first n-1 . having
c  w(1,n)+diag(n) .le. diag(n) signifies that the n-th function is lin-
c  early dependent to machine accuracy on the first n-1 functions, there
c  fore can safely be left out from the basis of approximating functions
c     the solution of a linear system
c                       c*x = b
c   is effected by the succession of the following  t w o  calls:
c     call bchfac ( w, nbands, nrow, diag )       , to get factorization
c     call bchslv ( w, nbands, nrow, b )          , to solve for x.
c
      integer nbands,nrow,   i,imax,j,jmax,n
      real w(nbands,nrow),diag(nrow),   ratio
      if (nrow .gt. 1)                  go to 9
      if (w(1,1) .gt. 0.) w(1,1) = 1./w(1,1)
                                        return
c                                        store diagonal of  c  in  diag.
    9 do 10 n=1,nrow
   10    diag(n) = w(1,n)
c                                                        factorization .
      do 20 n=1,nrow
         if (w(1,n)+diag(n) .gt. diag(n)) go to 15
         do 14 j=1,nbands
   14       w(j,n) = 0.
                                        go to 20
   15    w(1,n) = 1./w(1,n)
         imax = min0(nbands-1,nrow - n)
         if (imax .lt. 1)               go to 20
         jmax = imax
         do 18 i=1,imax
            ratio = w(i+1,n)*w(1,n)
            do 17 j=1,jmax
   17          w(j,n+i) = w(j,n+i) - w(j+i,n)*ratio
            jmax = jmax - 1
   18       w(i+1,n) = ratio
   20    continue
                                        return
      end