smooth Function

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.

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)).

16 February 2007

Carl DeBoor

Carl DeBoor,
A Practical Guide to Splines,
Springer, 2001,
ISBN: 0387953663,
LC: QA1.A647.v27.

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.