smooth Function

function smooth ( x, y, dy, npoint, s, v, a )

  !*****************************************************************************80
  !
  !! SMOOTH constructs the cubic smoothing spline to given data.
  !
  !  Discussion:
  !
  !    The data is of the form
  !
  !      ( X(1:NPOINT), Y(1:NPOINT) )
  !
  !    The cubic smoothing spline has as small a second derivative as
  !    possible, while
  !
  !      S(F) <= S,
  !
  !    where
  !
  !      S(F) = sum ( 1 <= I <= NPOINT ) ( ( ( Y(I) - F(X(I)) ) / DY(I) )**2.
  !
  !  Method:
  !
  !    The matrices Q' * D and Q' * D**2 * Q are constructed in SETUPQ from
  !    X and DY, as is the vector QTY = Q' * Y.
  !
  !    Then, for given P, the vector U is determined in CHOL1D as
  !    the solution of the linear system
  !
  !      ( 6 * (1-P) * Q' * D**2 * Q + P * R ) * U = QTY.
  !
  !    From U and this choice of smoothing parameter P, the smoothing spline F
  !    is obtained in the sense that:
  !
  !             F(X(.)) = Y - 6 (1-P) D**2 * Q * U,
  !      (d**2) F(X(.)) = 6 * P * U.
  !
  !    The smoothing parameter P is found, if possible, so that
  !
  !      SF(P) = S,
  !
  !    with SF(P) = S(F), where F is the smoothing spline as it depends
  !    on P.  If S = 0, then P = 1.  If SF(0) <= S, then P = 0.
  !    Otherwise, the secant method is used to locate an appropriate P in
  !    the open interval (0,1).
  !
  !    Specifically,
  !
  !      P(0) = 0,  P(1) = ( S - SF(0) ) / DSF
  !
  !    with
  !
  !      DSF = -24 * U' * R * U
  !
  !    a good approximation to
  !
  !      D(SF(0)) = DSF + 60 * (D*Q*U)' * (D*Q*U),
  !
  !    and U as obtained for P = 0.
  !
  !    After that, for N = 1, 2,...  until SF(P(N)) <= 1.01 * S, do:
  !    determine P(N+1) as the point at which the secant to SF at the
  !    points P(N) and P(N-1) takes on the value S.
  !
  !    If 1 <= P(N+1), choose instead P(N+1) as the point at which
  !    the parabola SF(P(N))*((1-.)/(1-P(N)))**2 takes on the value S.
  !
  !    Note that, in exact arithmetic, it is always the case that
  !      P(N+1) < P(N),
  !    hence
  !      SF(P(N+1)) < SF(P(N)).
  !
  !    Therefore, also stop the iteration, with final P = 1, in case
  !      SF(P(N)) <= SF(P(N+1)).
  !
  !  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 ) X(NPOINT), the abscissas, assumed to be
  !    strictly increasing.
  !
  !    Input, real ( kind = 8 ) Y(NPOINT), the corresponding ordinates.
  !
  !    Input, real ( kind = 8 ) DY(NPOINT), the data uncertainty estimates,
  !    which are assumed to be positive.
  !
  !    Input, integer ( kind = 4 ) NPOINT, the number of data points.
  !
  !    Input, real ( kind = 8 ) S, an upper bound on the discrete weighted mean
  !    square distance of the approximation F from the data.
  !
  !    Workspace, real ( kind = 8 ) V(NPOINT,7).
  !
  !    Workspace, real ( kind = 8 ) A(NPOINT,4).
  !
  !    Output, real ( kind = 8 ) A(NPOINT,4).
  !    A(*,1).....contains the sequence of smoothed ordinates.
  !    A(I,J) = d**(J-1) F(X(I)), for J = 2:4, I = 1:NPOINT-1.
  !    That is, the first three derivatives of the smoothing spline F at the
  !    left end of each of the data intervals.  Note that A would have to
  !    be transposed before it could be used in PPVALU.
  !
  !    Output, real ( kind = 8 ) SMOOTH, the value of the smoothing parameter.
  !
  implicit none

  integer ( kind = 4 ) npoint

  real ( kind = 8 ) a(npoint,4)
  real ( kind = 8 ) change
  real ( kind = 8 ) dy(npoint)
  integer ( kind = 4 ) i
  real ( kind = 8 ) oosf
  real ( kind = 8 ) ooss
  real ( kind = 8 ) p
  real ( kind = 8 ) prevq
  real ( kind = 8 ) prevsf
  real ( kind = 8 ) q
  real ( kind = 8 ) s
  real ( kind = 8 ) sfq
  real ( kind = 8 ) smooth
  real ( kind = 8 ) utru
  real ( kind = 8 ) v(npoint,7)
  real ( kind = 8 ) x(npoint)
  real ( kind = 8 ) y(npoint)

  call setupq ( x, dy, y, npoint, v, a(1,4) )

  if ( s <= 0.0D+00 ) then

     p = 1.0D+00
     call chol1d ( p, v, a(1,4), npoint, 1, a(1,3), a(1,1) )
     sfq = 0.0D+00

  else

     p = 0.0D+00
     call chol1d ( p, v, a(1,4), npoint, 1, a(1,3), a(1,1) )

     sfq = 36.0D+00 * dot_product ( a(1:npoint,1)**2, dy(1:npoint)**2 )

     if ( s < sfq ) then

        utru = 0.0D+00
        do i = 2, npoint
           utru = utru + v(i-1,4) * ( a(i-1,3) * ( a(i-1,3) + a(i,3) ) &
                + a(i,3)**2 )
        end do

        ooss = 1.0D+00 / sqrt ( s )
        oosf = 1.0D+00 / sqrt ( sfq )
        q = - ( oosf - ooss ) * sfq / ( 6.0D+00 * utru * oosf )
        !
        !  Secant iteration for the determination of P starts here.
        !
        prevq = 0.0D+00
        prevsf = oosf

        do

           call chol1d ( q / ( 1.0D+00 + q ), v, a(1,4), npoint, 1, &
                a(1,3), a(1,1) )

           sfq = 36.0D+00 * dot_product ( a(1:npoint,1)**2, dy(1:npoint)**2 ) &
                / ( 1.0D+00 + q )**2

           if ( abs ( sfq - s ) <= 0.01D+00 * s ) then
              exit
           end if

           oosf = 1.0D+00 / sqrt ( sfq )
           change = ( q - prevq ) / ( oosf - prevsf ) * ( oosf - ooss )
           prevq = q
           q = q - change
           prevsf = oosf

        end do

        p = q / ( 1.0D+00 + q )

     end if

  end if
  !
  !  Correct value of P has been found.
  !  Compute polynomial coefficients from Q * U in A(.,1).
  !
  smooth = sfq

  a(1:npoint,1) = y(1:npoint) - 6.0D+00 * ( 1.0D+00 - p ) &
       * dy(1:npoint)**2 * a(1:npoint,1)

  a(1:npoint,3) = a(1:npoint,3) * 6.0D+00 * p

  do i = 1, npoint - 1
     a(i,4) = ( a(i+1,3) - a(i,3) ) / v(i,4)
     a(i,2) = ( a(i+1,1) - a(i,1) ) / v(i,4) &
          - ( a(i,3) + a(i,4) / 3.0D+00 * v(i,4) ) / 2.0D+00 * v(i,4)
  end do

  return
end function smooth