ADAS Subroutine xxnbaf
SUBROUTINE XXNBAF (M, NCAP7 , X, Y, W, LAMDA , B, A, DIAG, C ,
& SS , IFAIL )
C
C
C-----------------------------------------------------------------------
C
C PURPOSE: Determines a least-square cubic spline approximation s(x)
C to the set of data points (x_r, y_r) with weights w_r,
C for r=1,2,...,m.
C
C The value of NCAP7 = ncap+7, where ncap is the number of
C intervals of the spline (number of interior knots + 1),
C and the values of the knots LAMDA(5), LAMDA(6), ...,
C LAMDA(NCAP7-4), interior to the data interval, are
C prescribed by the user.
C
C s has the property that it minimises ss, the sum of the
C squares of the weighted residuals eps(r)
C
C eps(r) = w(r) * ( s(x(r)) - y(r) ).
C
C The procedure produces the minimising value of ss and
C the coefficients c(1), c(2),...,c(q), where q=ncap+3=NCAP7-4,
C in the B-spline representation
C
C s(x) = c(1)*N1(x) + c(2)*N2(x) + ... + c(q)*Nq(x) .
C
C Here Ni(x) (i=1,2,...,q) denotes the normalised B-spline
C of degree 3 defined upon the knots lamda(i-4), lamda(i-3),
C lamda(i-2), lamda(i-1), and lamda(i).
C
C CALLING PROGRAM: VARIOUS
C
C SUBROUTINE:
C
C INPUT: (I*4) M = The number of data points.
C CONSTRAINT: M >= MDIST >= 4, where
C MDIST is the number of distinct x
C values in the data.
C
C INPUT: (I*4) NCAP7 = NBAR+7, where NBAR is the number of
C intervals of the spline (number
C of interior knots +1, i.e. the knots
C strictly in the range X(1) to X(M))
C over which the spline is defined.
C CONSTRAINT: 8<= NCAP7 <= MDIST+4,
C where MDIST is the number of distinct
C x values in the data.
C
C INPUT: (R*8) X() = The values x_r of the independent variable
C (abscissa), for r=1,2,...,m.
C CONSTRAINT:
C X(1)<=X(2)<=...<=X(M)
C
C INPUT: (R*8) Y() = The values y_r of the dependent variable
C (ordinate), for r=1,2,...,m.
C
C INPUT: (R*8) W() = The values w_r of the weights,
C for r=1,2,...,m.
C
C INPUT: (R*8) LAMDA() = LAMDA(i) must be set to the (i-4)th
C (interior) knot, i=5,6,...,nbar+3.
C CONSTRAINT:
C X(1) < LAMDA(5) <= LAMDA(6) ... <=
C ... <=LAMDA(NCAP7-4) < X(M) .
C
C INPUT: (I*8) IFAIL = 0 : stop if any error
C = 1 : continue if non-fatal error.
C
C OUTPUT: (R*8) LAMDA() = Input values are unchanged, and
C LAMDA(i), for i=1,2,3,4,NCAP7-3,
C NCAP7-2,NCAP7-1,NCAP7 contains the
C additional exterior knots introduced by
C the routine.
C
C
C OUTPUT: (R*8) C() = The coefficients of the B-spline N_i(x),
C for i=1,2,...,nbar+3. The remaining
C elements (from NBAR+4 to NBAR+7) are not
C used.
C
C OUTPUT: (R*8) SS = The residual sum of sqaures
C
C OUTPUT: (I*4) IFAIL = 0 : no error detected
C = 1 : the knots fail to satisfy the condition
C X(1) < LAMDA(5) <= LAMDA(6) <=...
C <= LAMDA(NCAP7-4) < X(M)
C = 2 : The weights are not strictly positive
C = 3 : The values of X(R), R=1,M are not in
C non-decreasing order.
C = 4 : NCAP7 < 8 (so that the number of
C interior knots is negative) or
C NCAP7 > MDIST + 4, where MDIST is the
C number of distinct x values in the data
C (so that there cannot be unique solution).
C = 5 : The conditions specified by Schoenberg
C and Whitney fail to hold for at least
C one subset of the distinct data abscissae.
C That is, there is no subset of NCAP7-4
C strictly increasing values,
C X(R(1)), X(R(2)),..., X(R(NCAP7-4)),
C among theabscissae such that
C
C X(R(1)) < LAMDA(1) < X(R(5))
C X(R(2)) < LAMDA(2) < X(R(6))
C ...
C X(R(NCAP7-8))<LAMDA(NCAP7-8)<X(R(NCAP7-4)).
C
C This means that there is no unique
C solution: there are regions containing
C too many knots compared with the
C number of data points.
C
C (R*8) B() = Set of distinct data abscissae
C
C (R*8) WORK2() = WORKSPACE
C
C (I*4) J = GENERAL INDEX
C (I*4) I = GENERAL INDEX
C (I*4) R = GENERAL INDEX
C (I*4) II = GENERAL INDEX
C (R*8) BI = GENERAL REAL
C (R*8) XI = GENERAL REAL
C
C ROUTINES: NONE
C
C AUTHORS: Alessandro C. Lanzafame, University of Strathclyde
C
C REFERENCE: Cox, M.G. and Hayes, J.G. "Curve fitting: A Guide and
C Suite of Algorithms for the Non-specialist User."
C Report NAC26, National Physical Laboratory, Middlessex,
C 1973.
C
C DATE: 12 January 1995
C
C VERSIION: 1.0a
C Alessandro Lanzafame, 12 January 1995.
C Directly derived from Algol text.
C (Error in passing woking variables)
C
C VERSION 1.0b
C Alessandro Lanzafame, 15 January 1995.
C DIAG(1:NCAP7-4) absorbed in matrix A(1:NCAP7-4,1).
C Matrix A(1:NCAP7-4,2:4) becomes A(1:NCAP7-4,1:4). This is to easy
C the passing of workspaces.
C WORK1 identified with B. WORK2 with A.
C Corrected index error in checking remaining Schoennber-Whitney
C conditions.
C COMPILING BUT NOT WORKING.
C
C VERSION 1.0c
C Alessandro Lanzafame, 15 January 1995.
C Knots shifted as in original Algol program text.
C
C UNIX-IDL PORT:
C
C VERSION: 1.1 DATE: 22-1-96
C MODIFIED: TIM HAMMOND (TESSELLA SUPPORT SERVICES PLC)
C - PUT UNDER SCCS CONTROL
C
C VERSION: 1.2 DATE: 06-07-2004
C MODIFIED: Allan Whiteford
C - Changed name from dxnbaf to xxnbaf.
C
C VERSION : 1.3 DATE: 10-04-2007
C MODIFIED : Allan Whiteford
C - Modified documentation as part of automated
C subroutine documentation preparation.
C
C-----------------------------------------------------------------------
C
C-----------------------------------------------------------------------
INTEGER IFAIL, M, NCAP7
REAL*8 A(1:NCAP7-4,2:4), B(M)
REAL*8 C(NCAP7), DIAG(1:NCAP7-4)
REAL*8 LAMDA(-3:NCAP7-4), SS, W(M)
REAL*8 X(M), Y(M)