splopt Subroutine

The optimal recovery scheme is of order K for data at TAU(1:N).

The interior knots T(K+1:N) are determined by Newton's method in 
such a way that the signum function which changes sign at
  T(K+1:N)  and nowhere else in ( TAU(1), TAU(N) ) is 
orthogonal to the spline space SPLINE ( K, TAU ) on that interval.

Let XI(J) be the current guess for T(K+J), J=1,...,N-K.  Then
the next Newton iterate is of the form

  XI(J) + (-)**(N-K-J)*X(J),  J=1,...,N-K,

with X the solution of the linear system

  C * X = D.

Here, for all J,

  C(I,J) = B(I)(XI(J)),

with B(I) the I-th B-spline of order K for the knot sequence TAU, 
for all I, and D is the vector given, for each I, by

  D(I) = sum ( -A(J), J=I,...,N ) * ( TAU(I+K) - TAU(I) ) / K,

with, for I = 1 to N-1:

  A(I) = sum ( (-)**(N-K-J)*B(I,K+1,TAU)(XI(J)), J=1,...,N-K )

and

  A(N) = -0.5.

See Chapter XIII of text and references there for a derivation.

The first guess for T(K+J) is sum ( TAU(J+1:J+K-1) ) / ( K - 1 ).

The iteration terminates if max ( abs ( X(J) ) ) < TOL, with

  TOL = TOLRTE * ( TAU(N) - TAU(1) ) / ( N - K ),

or else after NEWTMX iterations.

14 February 2007

Carl DeBoor

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

Input, real ( kind = 8 ) TAU(N), the interpolation points.
assumed to be nondecreasing, with TAU(I) < TAU(I+K), for all I.

Input, integer ( kind = 4 ) N, the number of data points.

Input, integer ( kind = 4 ) K, the order of the optimal recovery scheme 
to be used.

Workspace, real ( kind = 8 ) SCRTCH((N-K)*(2*K+3)+5*K+3).  The various
contents are specified in the text below.

Output, real ( kind = 8 ) T(N+K), the optimal knots ready for
use in optimal recovery.  Specifically, T(1:K) = TAU(1), 
T(N+1:N+K) = TAU(N), while the N - K interior knots T(K+1:N) 
are calculated.

Output, integer ( kind = 4 ) IFLAG, error indicator.
= 1, success.  T contains the optimal knots.
= 2, failure.  K < 3 or N < K or the linear system was singular.