bchfac Subroutine

subroutine bchfac ( w, nbands, nrow, diag )

  !*****************************************************************************80
  !
  !! BCHFAC constructs a Cholesky factorization of a matrix.
  !
  !  Discussion:
  !
  !    The factorization has the form
  !
  !      C = L * D * L'
  !  
  !    with L unit lower triangular and D diagonal, for a given matrix C of 
  !    order NROW, where C is symmetric positive semidefinite and banded, 
  !    having NBANDS diagonals at and below the main diagonal.
  ! 
  !    Gauss elimination is used, adapted to the symmetry and bandedness of C.
  ! 
  !    Near-zero pivots are handled in a special way.  The diagonal 
  !    element C(N,N) = W(1,N) is saved initially in DIAG(N), all N. 
  ! 
  !    At the N-th elimination step, the current pivot element, W(1,N), 
  !    is compared with its original value, DIAG(N).  If, as the result 
  !    of prior elimination steps, this element has been reduced by about 
  !    a word length, that is, if W(1,N) + DIAG(N) <= DIAG(N), then the pivot 
  !    is declared to be zero, and the entire N-th row is declared to
  !    be linearly dependent on the preceding rows.  This has the effect 
  !    of producing X(N) = 0 when solving C * X = B for X, regardless of B.
  ! 
  !    Justification for this is as follows.  In contemplated applications 
  !    of this program, the given equations are the normal equations for 
  !    some least-squares approximation problem, DIAG(N) = C(N,N) gives 
  !    the norm-square of the N-th basis function, and, at this point, 
  !    W(1,N) contains the norm-square of the error in the least-squares 
  !    approximation to the N-th basis function by linear combinations 
  !    of the first N-1.  
  !
  !    Having W(1,N)+DIAG(N) <= DIAG(N) signifies that the N-th function 
  !    is linearly dependent to machine accuracy on the first N-1 
  !    functions, therefore can safely be left out from the basis of 
  !    approximating functions.
  !
  !    The solution of a linear system C * X = B is effected by the 
  !    succession of the following two calls:
  ! 
  !      call bchfac ( w, nbands, nrow, diag )
  !
  !      call bchslv ( w, nbands, nrow, b, x )
  !
  !  Modified:
  !
  !    14 February 2007
  !
  !  Author:
  !
  !    Carl DeBoor
  !
  !  Reference:
  !
  !    Carl DeBoor,
  !    A Practical Guide to Splines,
  !    Springer, 2001,
  !    ISBN: 0387953663,
  !    LC: QA1.A647.v27.
  !
  !  Parameters:
  !
  !    Input/output, real ( kind = 8 ) W(NBANDS,NROW).
  !    On input, W contains the NBANDS diagonals in its rows, 
  !    with the main diagonal in row 1.  Precisely, W(I,J) 
  !    contains C(I+J-1,J), I=1,...,NBANDS, J=1,...,NROW.
  !    For example, the interesting entries of a seven diagonal
  !    symmetric matrix C of order 9 would be stored in W as
  ! 
  !      11 22 33 44 55 66 77 88 99
  !      21 32 43 54 65 76 87 98  *
  !      31 42 53 64 75 86 97  *  *
  !      41 52 63 74 85 96  *  *  *
  !
  !    Entries of the array not associated with an
  !    entry of C are never referenced.
  !    On output, W contains the Cholesky factorization 
  !    C = L*D*L', with W(1,I) containing 1/D(I,I) and W(I,J) 
  !    containing L(I-1+J,J), I=2,...,NBANDS.
  !
  !    Input, integer ( kind = 4 ) NBANDS, indicates the bandwidth of the
  !    matrix C, that is, C(I,J) = 0 for NBANDS < abs(I-J).
  ! 
  !    Input, integer ( kind = 4 ) NROW, is the order of the matrix C.
  ! 
  !    Work array, real ( kind = 8 ) DIAG(NROW).
  !
  implicit none

  integer ( kind = 4 ) nbands
  integer ( kind = 4 ) nrow

  real ( kind = 8 ) diag(nrow)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) imax
  integer ( kind = 4 ) j
  integer ( kind = 4 ) jmax
  integer ( kind = 4 ) n
  real ( kind = 8 ) ratio
  real ( kind = 8 ) w(nbands,nrow)

  if ( nrow <= 1 ) then
     if ( 0.0D+00 < w(1,1) ) then
        w(1,1) = 1.0D+00 / w(1,1)
     end if
     return
  end if
  !
  !  Store the diagonal.
  !
  diag(1:nrow) = w(1,1:nrow)
  !
  !  Factorization.
  !
  do n = 1, nrow

     if ( w(1,n) + diag(n) <= diag(n) ) then
        w(1:nbands,n) = 0.0D+00
     else

        w(1,n) = 1.0D+00 / w(1,n)

        imax = min ( nbands - 1, nrow - n )

        jmax = imax

        do i = 1, imax

           ratio = w(i+1,n) * w(1,n)

           do j = 1, jmax
              w(j,n+i) = w(j,n+i) - w(j+i,n) * ratio
           end do

           jmax = jmax-1
           w(i+1,n) = ratio

        end do

     end if

  end do

  return
end subroutine bchfac