chol1d Subroutine

subroutine chol1d ( p, v, qty, npoint, ncol, u, qu )

  !*****************************************************************************80
  !
  !! CHOL1D sets up and solves linear systems needed by SMOOTH.
  !
  !  Discussion:
  !
  !    This routine constructs the upper three diagonals of
  !
  !      V(I,J), I = 2 to NPOINT-1, J=1,3, 
  !
  !    of the matrix 
  !
  !      6 * (1-P) * Q' * (D**2) * Q + P * R.
  ! 
  !    It then computes its L*L' decomposition and stores it also
  !    in V, then applies forward and back substitution to the right hand side 
  !
  !      Q'*Y 
  !
  !    in QTY to obtain the solution in U.
  !
  !  Modified:
  !
  !    16 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 ) P, the smoothing parameter that defines 
  !    the linear system.
  !
  !    Input/output, real ( kind = 8 ) V(NPOINT,7), contains data used
  !    to define the linear system, some of which is determined by
  !    routine SETUPQ.
  !
  !    Input, real ( kind = 8 ) QTY(NPOINT), the value of Q' * Y.
  !
  !    Input, integer ( kind = 4 ) NPOINT, the number of equations.
  !
  !    Input, integer ( kind = 4 ) NCOL, an unused parameter, which may be 
  !    set to 1.
  !
  !    Output, real ( kind = 8 ) U(NPOINT), the solution.
  !
  !    Output, real ( kind = 8 ) QU(NPOINT), the value of Q * U.
  !
  implicit none

  integer ( kind = 4 ) npoint

  integer ( kind = 4 ) i
  integer ( kind = 4 ) ncol
  real ( kind = 8 ) p
  real ( kind = 8 ) qty(npoint)
  real ( kind = 8 ) qu(npoint)
  real ( kind = 8 ) u(npoint)
  real ( kind = 8 ) v(npoint,7)
  real ( kind = 8 ) prev
  real ( kind = 8 ) ratio
  real ( kind = 8 ) six1mp
  real ( kind = 8 ) twop
  !
  !  Construct 6*(1-P)*Q'*(D**2)*Q + P*R.
  !
  six1mp = 6.0D+00 * ( 1.0D+00 - p )
  twop = 2.0D+00 * p

  v(2:npoint-1,1) = six1mp * v(2:npoint-1,5) &
       + twop * ( v(1:npoint-2,4) + v(2:npoint-1,4) )
  v(2:npoint-1,2) = six1mp * v(2:npoint-1,6) + p * v(2:npoint-1,4)
  v(2:npoint-1,3) = six1mp * v(2:npoint-1,7)

  if ( npoint < 4 ) then

     u(1) = 0.0D+00
     u(2) = qty(2) / v(2,1)
     u(3) = 0.0D+00
     !
     !  Factorization.
     !
  else

     do i = 2, npoint-2
        ratio = v(i,2) / v(i,1)
        v(i+1,1) = v(i+1,1) - ratio * v(i,2)
        v(i+1,2) = v(i+1,2) - ratio * v(i,3)
        v(i,2) = ratio
        ratio = v(i,3) / v(i,1)
        v(i+2,1) = v(i+2,1) - ratio * v(i,3)
        v(i,3) = ratio
     end do
     !
     !  Forward substitution
     !
     u(1) = 0.0D+00
     v(1,3) = 0.0D+00
     u(2) = qty(2)
     do i = 2, npoint-2
        u(i+1) = qty(i+1) - v(i,2) * u(i) - v(i-1,3) * u(i-1)
     end do
     !
     !  Back substitution.
     !
     u(npoint) = 0.0D+00
     u(npoint-1) = u(npoint-1) / v(npoint-1,1)

     do i = npoint-2, 2, -1
        u(i) = u(i) / v(i,1) - u(i+1) * v(i,2) - u(i+2) * v(i,3)
     end do

  end if
  !
  !  Construct Q * U.
  !
  prev = 0.0D+00
  do i = 2, npoint
     qu(i) = ( u(i) - u(i-1) ) / v(i-1,4)
     qu(i-1) = qu(i) - prev
     prev = qu(i)
  end do

  qu(npoint) = -qu(npoint)

  return
end subroutine chol1d