bchslv Subroutine

subroutine bchslv ( w, nbands, nrow, b )

  !*****************************************************************************80
  !
  !! BCHSLV solves a banded symmetric positive definite system.
  !
  !  Discussion:
  !
  !    The system is of the form:
  !
  !      C * X = B 
  !  
  !    and the Cholesky factorization of C has been constructed 
  !    by BCHFAC.
  ! 
  !    With the factorization 
  !
  !      C = L * D * L'
  !
  !    available, where L is unit lower triangular and D is diagonal, 
  !    the triangular system 
  !
  !      L * Y = B 
  !
  !    is solved for Y (forward substitution), Y is stored in B, the 
  !    vector D**(-1)*Y is computed and stored in B, then the 
  !    triangular system L'*X = D**(-1)*Y is solved for X 
  !    (back substitution).
  !
  !  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, real ( kind = 8 ) W(NBANDS,NROW), the Cholesky factorization for C, 
  !    as computed by BCHFAC.
  ! 
  !    Input, integer ( kind = 4 ) NBANDS, the bandwidth of C.
  !
  !    Input, integer ( kind = 4 ) NROW, the order of the matrix C.
  ! 
  !    Input/output, real ( kind = 8 ) B(NROW).
  !    On input, the right hand side.
  !    On output, the solution.
  !
  implicit none

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

  real ( kind = 8 ) b(nrow)
  integer ( kind = 4 ) j
  integer ( kind = 4 ) n
  real ( kind = 8 ) w(nbands,nrow)

  if ( nrow <= 1 ) then
     b(1) = b(1) * w(1,1)
     return
  end if
  !
  !  Forward substitution. 
  !  Solve L*Y = B.
  !
  do n = 1, nrow

     do j = 1, min ( nbands - 1, nrow - n )
        b(j+n) = b(j+n) - w(j+1,n) * b(n)
     end do

  end do
  !
  !  Back substitution. 
  !  Solve L'*X = D**(-1)*Y.
  !
  do n = nrow, 1, -1

     b(n) = b(n) * w(1,n)

     do j = 1, min ( nbands - 1, nrow - n )
        b(n) = b(n) - w(j+1,n) * b(j+n)
     end do

  end do

  return
end subroutine bchslv