1*DECK DEFC
2      SUBROUTINE DEFC (NDATA, XDATA, YDATA, SDDATA, NORD, NBKPT, BKPT,
3     +   MDEIN, MDEOUT, COEFF, LW, W)
4C***BEGIN PROLOGUE  DEFC
5C***PURPOSE  Fit a piecewise polynomial curve to discrete data.
6C            The piecewise polynomials are represented as B-splines.
7C            The fitting is done in a weighted least squares sense.
8C***LIBRARY   SLATEC
9C***CATEGORY  K1A1A1, K1A2A, L8A3
10C***TYPE      DOUBLE PRECISION (EFC-S, DEFC-D)
11C***KEYWORDS  B-SPLINE, CONSTRAINED LEAST SQUARES, CURVE FITTING
12C***AUTHOR  Hanson, R. J., (SNLA)
13C***DESCRIPTION
14C
15C      This subprogram fits a piecewise polynomial curve
16C      to discrete data.  The piecewise polynomials are
17C      represented as B-splines.
18C      The fitting is done in a weighted least squares sense.
19C
20C      The data can be processed in groups of modest size.
21C      The size of the group is chosen by the user.  This feature
22C      may be necessary for purposes of using constrained curve fitting
23C      with subprogram DFC( ) on a very large data set.
24C
25C      For a description of the B-splines and usage instructions to
26C      evaluate them, see
27C
28C      C. W. de Boor, Package for Calculating with B-Splines.
29C                     SIAM J. Numer. Anal., p. 441, (June, 1977).
30C
31C      For further discussion of (constrained) curve fitting using
32C      B-splines, see
33C
34C      R. J. Hanson, Constrained Least Squares Curve Fitting
35C                   to Discrete Data Using B-Splines, a User's
36C                   Guide. Sandia Labs. Tech. Rept. SAND-78-1291,
37C                   December, (1978).
38C
39C  Input.. All TYPE REAL variables are DOUBLE PRECISION
40C      NDATA,XDATA(*),
41C      YDATA(*),
42C      SDDATA(*)
43C                         The NDATA discrete (X,Y) pairs and the Y value
44C                         standard deviation or uncertainty, SD, are in
45C                         the respective arrays XDATA(*), YDATA(*), and
46C                         SDDATA(*).  No sorting of XDATA(*) is
47C                         required.  Any non-negative value of NDATA is
48C                         allowed.  A negative value of NDATA is an
49C                         error.  A zero value for any entry of
50C                         SDDATA(*) will weight that data point as 1.
51C                         Otherwise the weight of that data point is
52C                         the reciprocal of this entry.
53C
54C      NORD,NBKPT,
55C      BKPT(*)
56C                         The NBKPT knots of the B-spline of order NORD
57C                         are in the array BKPT(*).  Normally the
58C                         problem data interval will be included between
59C                         the limits BKPT(NORD) and BKPT(NBKPT-NORD+1).
60C                         The additional end knots BKPT(I),I=1,...,
61C                         NORD-1 and I=NBKPT-NORD+2,...,NBKPT, are
62C                         required to compute the functions used to fit
63C                         the data.  No sorting of BKPT(*) is required.
64C                         Internal to DEFC( ) the extreme end knots may
65C                         be reduced and increased respectively to
66C                         accommodate any data values that are exterior
67C                         to the given knot values.  The contents of
68C                         BKPT(*) is not changed.
69C
70C                         NORD must be in the range 1 .LE. NORD .LE. 20.
71C                         The value of NBKPT must satisfy the condition
72C                         NBKPT .GE. 2*NORD.
73C                         Other values are considered errors.
74C
75C                         (The order of the spline is one more than the
76C                         degree of the piecewise polynomial defined on
77C                         each interval.  This is consistent with the
78C                         B-spline package convention.  For example,
79C                         NORD=4 when we are using piecewise cubics.)
80C
81C        MDEIN
82C                         An integer flag, with one of two possible
83C                         values (1 or 2), that directs the subprogram
84C                         action with regard to new data points provided
85C                         by the user.
86C
87C                         =1  The first time that DEFC( ) has been
88C                         entered.  There are NDATA points to process.
89C
90C                         =2  This is another entry to DEFC().  The sub-
91C                         program DEFC( ) has been entered with MDEIN=1
92C                         exactly once before for this problem.  There
93C                         are NDATA new additional points to merge and
94C                         process with any previous points.
95C                         (When using DEFC( ) with MDEIN=2 it is import-
96C                         ant that the set of knots remain fixed at the
97C                         same values for all entries to DEFC( ).)
98C       LW
99C                         The amount of working storage actually
100C                         allocated for the working array W(*).
101C                         This quantity is compared with the
102C                         actual amount of storage needed in DEFC( ).
103C                         Insufficient storage allocated for W(*) is
104C                         an error.  This feature was included in DEFC
105C                         because misreading the storage formula
106C                         for W(*) might very well lead to subtle
107C                         and hard-to-find programming bugs.
108C
109C                         The length of the array W(*) must satisfy
110C
111C                         LW .GE. (NBKPT-NORD+3)*(NORD+1)+
112C                                 (NBKPT+1)*(NORD+1)+
113C                               2*MAX(NDATA,NBKPT)+NBKPT+NORD**2
114C
115C  Output.. All TYPE REAL variables are DOUBLE PRECISION
116C      MDEOUT
117C                         An output flag that indicates the status
118C                         of the curve fit.
119C
120C                         =-1  A usage error of DEFC( ) occurred.  The
121C                         offending condition is noted with the SLATEC
122C                         library error processor, XERMSG( ).  In case
123C                         the working array W(*) is not long enough, the
124C                         minimal acceptable length is printed.
125C
126C                         =1  The B-spline coefficients for the fitted
127C                         curve have been returned in array COEFF(*).
128C
129C                         =2  Not enough data has been processed to
130C                         determine the B-spline coefficients.
131C                         The user has one of two options.  Continue
132C                         to process more data until a unique set
133C                         of coefficients is obtained, or use the
134C                         subprogram DFC( ) to obtain a specific
135C                         set of coefficients.  The user should read
136C                         the usage instructions for DFC( ) for further
137C                         details if this second option is chosen.
138C      COEFF(*)
139C                         If the output value of MDEOUT=1, this array
140C                         contains the unknowns obtained from the least
141C                         squares fitting process.  These N=NBKPT-NORD
142C                         parameters are the B-spline coefficients.
143C                         For MDEOUT=2, not enough data was processed to
144C                         uniquely determine the B-spline coefficients.
145C                         In this case, and also when MDEOUT=-1, all
146C                         values of COEFF(*) are set to zero.
147C
148C                         If the user is not satisfied with the fitted
149C                         curve returned by DEFC( ), the constrained
150C                         least squares curve fitting subprogram DFC( )
151C                         may be required.  The work done within DEFC( )
152C                         to accumulate the data can be utilized by
153C                         the user, if so desired.  This involves
154C                         saving the first (NBKPT-NORD+3)*(NORD+1)
155C                         entries of W(*) and providing this data
156C                         to DFC( ) with the "old problem" designation.
157C                         The user should read the usage instructions
158C                         for subprogram DFC( ) for further details.
159C
160C  Working Array.. All TYPE REAL variables are DOUBLE PRECISION
161C      W(*)
162C                         This array is typed DOUBLE PRECISION.
163C                         Its length is  specified as an input parameter
164C                         in LW as noted above.  The contents of W(*)
165C                         must not be modified by the user between calls
166C                         to DEFC( ) with values of MDEIN=1,2,2,... .
167C                         The first (NBKPT-NORD+3)*(NORD+1) entries of
168C                         W(*) are acceptable as direct input to DFC( )
169C                         for an "old problem" only when MDEOUT=1 or 2.
170C
171C  Evaluating the
172C  Fitted Curve..
173C                         To evaluate derivative number IDER at XVAL,
174C                         use the function subprogram DBVALU( ).
175C
176C                         F = DBVALU(BKPT,COEFF,NBKPT-NORD,NORD,IDER,
177C                                      XVAL,INBV,WORKB)
178C
179C                         The output of this subprogram will not be
180C                         defined unless an output value of MDEOUT=1
181C                         was obtained from DEFC( ), XVAL is in the data
182C                         interval, and IDER is nonnegative and .LT.
183C                         NORD.
184C
185C                         The first time DBVALU( ) is called, INBV=1
186C                         must be specified.  This value of INBV is the
187C                         overwritten by DBVALU( ).  The array WORKB(*)
188C                         must be of length at least 3*NORD, and must
189C                         not be the same as the W(*) array used in the
190C                         call to DEFC( ).
191C
192C                         DBVALU( ) expects the breakpoint array BKPT(*)
193C                         to be sorted.
194C
195C***REFERENCES  R. J. Hanson, Constrained least squares curve fitting
196C                 to discrete data using B-splines, a users guide,
197C                 Report SAND78-1291, Sandia Laboratories, December
198C                 1978.
199C***ROUTINES CALLED  DEFCMN
200C***REVISION HISTORY  (YYMMDD)
201C   800801  DATE WRITTEN
202C   890531  Changed all specific intrinsics to generic.  (WRB)
203C   890531  REVISION DATE from Version 3.2
204C   891214  Prologue converted to Version 4.0 format.  (BAB)
205C   900510  Change Prologue comments to refer to XERMSG.  (RWC)
206C   900607  Editorial changes to Prologue to make Prologues for EFC,
207C           DEFC, FC, and DFC look as much the same as possible.  (RWC)
208C   920501  Reformatted the REFERENCES section.  (WRB)
209C***END PROLOGUE  DEFC
210C
211C      SUBROUTINE           FUNCTION/REMARKS
212C
213C      DBSPVN( )          COMPUTE FUNCTION VALUES OF B-SPLINES.  FROM
214C                         THE B-SPLINE PACKAGE OF DE BOOR NOTED ABOVE.
215C
216C      DBNDAC( ),         BANDED LEAST SQUARES MATRIX PROCESSORS.
217C      DBNDSL( )          FROM LAWSON-HANSON, SOLVING LEAST
218C                         SQUARES PROBLEMS.
219C
220C      DSORT( )           DATA SORTING SUBROUTINE, FROM THE
221C                         SANDIA MATH. LIBRARY, SAND77-1441.
222C
223C      XERMSG( )          ERROR HANDLING ROUTINE
224C                         FOR THE SLATEC MATH. LIBRARY.
225C                         SEE SAND78-1189, BY R. E. JONES.
226C
227C      DCOPY( ),DSCAL( )  SUBROUTINES FROM THE BLAS PACKAGE.
228C
229C                         WRITTEN BY R. HANSON, SANDIA NATL. LABS.,
230C                         ALB., N. M., AUGUST-SEPTEMBER, 1980.
231C
232      DOUBLE PRECISION BKPT(*),COEFF(*),W(*),SDDATA(*),XDATA(*),YDATA(*)
233      INTEGER LW, MDEIN, MDEOUT, NBKPT, NDATA, NORD
234C
235      EXTERNAL DEFCMN
236C
237      INTEGER LBF, LBKPT, LG, LPTEMP, LWW, LXTEMP, MDG, MDW
238C
239C***FIRST EXECUTABLE STATEMENT  DEFC
240C     LWW=1               USAGE IN DEFCMN( ) OF W(*)..
241C     LWW,...,LG-1        W(*,*)
242C
243C     LG,...,LXTEMP-1     G(*,*)
244C
245C     LXTEMP,...,LPTEMP-1 XTEMP(*)
246C
247C     LPTEMP,...,LBKPT-1  PTEMP(*)
248C
249C     LBKPT,...,LBF       BKPT(*) (LOCAL TO DEFCMN( ))
250C
251C     LBF,...,LBF+NORD**2 BF(*,*)
252C
253      MDG = NBKPT+1
254      MDW = NBKPT-NORD+3
255      LWW = 1
256      LG = LWW + MDW*(NORD+1)
257      LXTEMP = LG + MDG*(NORD+1)
258      LPTEMP = LXTEMP + MAX(NDATA,NBKPT)
259      LBKPT  = LPTEMP + MAX(NDATA,NBKPT)
260      LBF = LBKPT + NBKPT
261      CALL DEFCMN(NDATA,XDATA,YDATA,SDDATA,
262     1         NORD,NBKPT,BKPT,
263     2         MDEIN,MDEOUT,
264     3         COEFF,
265     4         W(LBF),W(LXTEMP),W(LPTEMP),W(LBKPT),
266     5         W(LG),MDG,W(LWW),MDW,LW)
267      RETURN
268      END
269