setupq Subroutine

subroutine setupq ( x, dx, y, npoint, v, qty )

  !*****************************************************************************80
  !
  !! SETUPQ is to be called in SMOOTH.
  !
  !  Discussion:
  !
  !    Put DELX = X(*+1) - X(*) into V(*,4).
  !
  !    Put the three bands of Q' * D into V(*,1:3).
  !
  !    Put the three bands of ( D * Q )' * ( D * Q ) at and above the diagonal
  !    into V(*,5:7).
  !
  !    Here, Q is the tridiagonal matrix of order ( NPOINT-2, NPOINT )
  !    with general row 
  !
  !      1/DELX(I), -1/DELX(I)-1/DELX(I+1), 1/DELX(I+1)
  !
  !    and D is the diagonal matrix with general row DX(I).
  !
  !  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 ) X(NPOINT), the abscissas, assumed to be 
  !    strictly increasing.
  !
  !    Input, real ( kind = 8 ) DX(NPOINT), the data uncertainty estimates,
  !    which are assumed to be positive.
  !
  !    Input, real ( kind = 8 ) Y(NPOINT), the corresponding ordinates.
  !
  !    Input, integer ( kind = 4 ) NPOINT, the number of data points.
  !
  !    Output, real ( kind = 8 ) V(NPOINT,7), contains data needed for
  !    the smoothing computation.
  !
  !    Output, real ( kind = 8 ) QTY(NPOINT), the value of Q' * Y.
  !
  implicit none

  integer ( kind = 4 ) npoint

  real ( kind = 8 ) diff
  real ( kind = 8 ) dx(npoint)
  integer ( kind = 4 ) i
  real ( kind = 8 ) prev
  real ( kind = 8 ) qty(npoint)
  real ( kind = 8 ) v(npoint,7)
  real ( kind = 8 ) x(npoint)
  real ( kind = 8 ) y(npoint)

  v(1:npoint-1,4) = x(2:npoint) - x(1:npoint-1)

  v(2:npoint-1,1) = dx(1:npoint-2) / v(1:npoint-2,4)  
  v(npoint,1) = 0.0D+00

  v(2:npoint-1,2) = - dx(2:npoint-1) / v(2:npoint-1,4) &
       - dx(2:npoint-1) / v(1:npoint-2,4)

  v(2:npoint-1,3) = dx(3:npoint) / v(2:npoint-1,4)

  v(2:npoint-1,5) = v(2:npoint-1,1)**2 &
       + v(2:npoint-1,2)**2 &
       + v(2:npoint-1,3)**2

  v(2:npoint-2,6) = v(2:npoint-2,2) * v(3:npoint-1,1) &
       + v(2:npoint-2,3) * v(3:npoint-1,2)
  v(npoint-1,6) = 0.0D+00

  v(2:npoint-3,7) = v(2:npoint-3,3) * v(4:npoint-1,1)   
  v(npoint-2,7) = 0.0D+00
  v(npoint-1,7) = 0.0D+00
  !
  !  Construct Q' * Y in QTY.
  !
  prev = ( y(2) - y(1) ) / v(1,4)
  do i = 2, npoint - 1
     diff = ( y(i+1) - y(i) ) / v(i,4)
     qty(i) = diff - prev
     prev = diff
  end do

  return
end subroutine setupq