/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:57 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include "dsfit.h"
#include <stdlib.h>
/* PARAMETER translations */
#define KMAX 20
#define ONE 1.0e0
#define ZERO 0.0e0
/* end of PARAMETER translations */
void /*FUNCTION*/ dsfit(
double xi[],
double yi[],
double sdi[],
long nxy,
long korder,
long nc,
double tknots[],
double bcoef[],
double *sigfac,
long *ierr1,
long ldw,
double *w)
{
#define W(I_,J_) (*(w+(I_)*(ldw)+(J_)))
LOGICAL32 direct, usewt1;
long int i, ierr2, ierr3, irnow, iseg, isgmax, j, jpoint, jpt,
jtprev, k, kordp1, ksize, left, nband, nt, nxyp1;
double dof, p[KMAX], rnorm, sdijp, wt, wt1;
/* OFFSET Vectors w/subscript range: 1 to dimension */
double *const Bcoef = &bcoef[0] - 1;
double *const P = &p[0] - 1;
double *const Sdi = &sdi[0] - 1;
double *const Tknots = &tknots[0] - 1;
double *const Xi = &xi[0] - 1;
double *const Yi = &yi[0] - 1;
/* end of OFFSET VECTORS */
/* Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
* ALL RIGHTS RESERVED.
* Based on Government Sponsored Research NAS7-03001.
*>> 2000-12-03 DSFIT Krogh Declared SDI(*) so ref. with SD(1) o.k.
*>> 1995-11-21 DSFIT Krogh Converted from SFTRAN to Fortran 77.
*>> 1994-10-19 DSFIT Krogh Changes to use M77CON
*>> 1994-01-31 DSFIT CLL Add test for SDI(i) .le. 0 when SDI(1) > 0.
*>> 1993-01-12 CLL Using DSBASD in place of _SBAS.
*>> 1992-12-15 CLL. Change error message No. 600
*>> 1992-11-18 CLL. Change order of last four arguments.
*>> 1992-11-17 CLL. Permit XI() nondecreasing or nonincreasing.
*>> 1992-10-27 C. L. Lawson, JPL
*>> 1988-03-21 C. L. Lawson, JPL
* Least squares fit to discrete data by a spline function of
* order KORDER. Note that the order is one greater than the degree
* of the polynomial pieces. Example: The order of a cubic spline
* is 4.
* TKNOTS() must contain NT values, called knots, indexed from
* 1 to NT, with NT = NC+KORDER. These values must be nondecreasing.
* Repeated values are permitted. The first and last KORDER-1 values
* in TKNOTS() are needed to support the deBoor method of
* representing splines. The "proper fitting interval" is from
* A = TKNOTS(KORDER) to B = TKNOTS(NT+1-KORDER). One acceptable way
* to set the first and last KORDER-1 knots is to set the first
* KORDER-1 to A and the last KORDER-1 to B.
* Continuity of the spline at knots interior to (A, B) will be of
* order KORDER-2, unless a knot is repeated, in which case the order
* of continuity will be decreased at that knot.
* ------------------------------------------------------------------
* The linear algebra methods were designed by C.L.Lawson and
* R.J.Hanson. The method of representing spline functions is due
* to Carl deBoor. References:
* "SOLVING LEAST SQUARES PROBLEMS", by Lawson and Hanson,
* Prentice-Hall, 1974.
* "A PRACTICAL GUIDE TO SPLINES" by Carl de Boor,
* Springer-Verlag, 1978.
* Programming and later changes and corrections by Lawson,Hanson,
* T.Lang, and D.Campbell, Sept 1968, Nov 1969, and Aug 1970.
* 1974 5/21, C.L.Lawson, Fixed bug that caused ISEG to get too big.
* Also changed to exit immidiately if TKNOTS() array is not
* strictly increasing.
* 1984 July 10. Modified for improved portability and to conform
* to Fortran 77. C. L. Lawson, JPL.
* Added calls to the error message subrs.
* 7/23/87 CLL. Added the IERR1 argument.
* March 1988, CLL. Introduced the use of deBoor's spline methods.
* 1992-11-17 CLL Permit XI() to be either nondecreasing or
* nonincreasing.
* ------------------------------------------------------------------
* SUBROUTINE ARGUMENTS
*
* (XI(i),i=1,NXY) [in] Abcissas of data to be fitted. Require
* this data be sorted: either nondecreasing or nonincreasing.
*
* (YI(i),i=1,NXY) [in] Ordinates of data to be fitted.
*
* (SDI(i),i=1,NXY) [in] User may use this array to assign an
* a priori standard deviation of error to each
* YI(i) value. The weighted fitting algorithm will take
* account of these. Optionally the user may set SDI(1) to
* a negative value. Then this subr will use ABS(SDI(1)) as
* the standard deviation for each YI(i) value. In this case
* the SDI() array can be dimensioned SDI(1).
* If SDI(1) = 0., the subr issues an error message and returns.
*
* NXY [in] No. of data pairs, (XI(i), YI(i)), and no. of elts
* in SDI() if SDI(1) is positive. Require NXY .ge. 4.
*
* KORDER [in] Order of the spline basis functions. Note that the
* polynomial degree of the spline segments is one less than the
* order. Example: the order of a cubic spline is 4.
* Require KORDER .ge. 1. Internal arrays in subroutines used put
* an upper limit of 20 on KORDER.
*
* NC [in] No. of coefficients to be determined. This is the
* number of degrees of freedom for the least squares problem.
* To have a nonsingular problem one must have
* NXY .ge. NC, and the distribution of the interior knots
* must not be too skewed relative to the data abcissas.
* If singularity is detected, an error message will be
* issued by the subr that solves the band matrix.
*
* (TKNOTS(j),j=1,NT, where NT = NC+KORDER) [in] This is the deBoor
* knot sequence for definition of the spline basis functions.
* See remarks above.
*
* BCOEF() [out] An array of length NC into which the computed
* coefficients defining the fitted curve will be stored. These
* are coeffients relative to B-spline basis functions. BCOEF(I)
* is associated with the basis function whose support interval
* runs from TKNOTS(I) to TKNOTS(I+KORDER).
*
* SIGFAC [out] The subr sets SIGFAC to RNORM / sqrt(DOF) where
* RNORM = sqrt( sum over i of [( (yfit(i) - YI(i))/SDI(i))**2])
* and DOF = max(1, NXY - NC)
*
* IERR1 [out] Error status indicator. Note that IERR2 comes from
* DBACC and IERR3 comes from DBSOL.
*
* = 0 means no errors detected. */
/* = 100 means NC .lt. 1 .or. NC .gt. NXY
* = 200 means TKNOTS(I) .gt. TKNOTS(I+1)
* = 250 means TKNOTS(I) .ge. TKNOTS(I+KORDER)
* = 300 means LDW .lt. NC+2
* = 400 means The XI's are not sorted.
* = 600 means Need larger dimension LDW.
* = 700 + IERR2 means IERR2 .ne. 0
* = 800 + IERR2 means IERR2 .ne. 0
* = 900 + IERR2 means IERR2 .ne. 0
* = 1000 + IERR3 means IERR3 .ne. 0 due to singularity
* detected in _BSOL.
* = 1100 means SDI(1) = zero.
* = 1200 means SDI(1) > zero and SDI(i) .le. zero for some i.
*
* LDW [in] Leading dimension of W(). Must satisfy LDW .ge. NC + 2
* Let IXMAX denote the max no. of data abcissas, XI(i),
* in any one knot interval, i.e. between TKNOTS(j) and
* TKNOTS(j+1) for some j. The subr will be more efficient
* if LDW is at least NC + 1 + IXMAX.
*
* W() [scratch] Work space dimensioned W(LDW,KORDER+1).
* ------------------------------------------------------------------
* Important internal variables.
*
* ISEG Index of current spline segment. ISEG runs from 1 to
* NC+1-KORDER. The knot interval associated with index ISEG
* is T(ISEG+KORDER-1). Note that the union of these
* segments is the "proper fitting interval".
* ISEG also tells the band matrix subroutine the column index
* of the least squares matrix with which the first col
* of the new block of data in G() is to be associated.
* There will be KORDER basis functions that are nonzero on
* segment ISEG. They are indexed from ISEG to ISEG+KORDER-1.
* KORDP1 = KORDER+1
* KSIZE Number of rows in current block.
* JPOINT Current data pointer.
* ------------------------------------------------------------------
*--D replaces "?": ?SFIT, ?BACC, ?BSOL, ?SBASD, ?ERV1
* Both versions use ERMSG, IERM1, IERV1
* Other subrs needed: DHTCC, ERMSG, ERFIN
* ------------------------------------------------------------------ */
/* ------------------------------------------------------------------ */
nband = korder;
kordp1 = korder + 1;
nt = nc + korder;
isgmax = nc + 1 - korder;
*ierr1 = 0;
/* Exit immediately if NC .lt. 1 or NXY .lt. NC or if the
* knots fail to be nondecreasing.
* */
if (nc < 1 || nxy < nc)
{
*ierr1 = 100;
ierm1( "DSFIT", *ierr1, 0, "Require NC .ge. 1 and NXY .ge. NC"
, "NC", nc, ',' );
ierv1( "NXY", nxy, '.' );
goto L_200;
}
if (korder > KMAX)
{
*ierr1 = 150;
ierm1( "DSFIT", *ierr1, 0, "Require KORDER .le. KMAX.", "KORDER"
, korder, ',' );
ierv1( "KMAX", KMAX, '.' );
goto L_200;
}
for (i = 1; i <= (nt - 1); i++)
{
if (Tknots[i] > Tknots[i + 1])
{
*ierr1 = 200;
ierm1( "DSFIT", *ierr1, 0, "Require knots, TKNOTS(I), to be nondecreasing."
, "I", i, ',' );
derv1( "TKNOTS(I)", Tknots[i], ',' );
derv1( "TKNOTS(I+1)", Tknots[i + 1], '.' );
goto L_200;
}
}
for (i = 1; i <= nc; i++)
{
if (Tknots[i] >= Tknots[i + korder])
{
*ierr1 = 250;
ierm1( "DSFIT", *ierr1, 0, "Require TKNOTS(I) < TKNOTS(I+KORDER)."
, "I", i, ',' );
derv1( "TKNOTS(I)", Tknots[i], ',' );
derv1( "TKNOTS(I+KORDER)", Tknots[i + korder], '.' );
goto L_200;
}
}
/* Require LDW .ge. NC+2
* */
if (ldw < nc + 2)
{
*ierr1 = 300;
ierm1( "DSFIT", *ierr1, 0, "Require LDW .ge. NC+2", "LDW",
ldw, ',' );
ierv1( "NC", nc, '.' );
goto L_200;
}
/* ------------------------------------------------------------------
* TEST SDI(1) */
if (Sdi[1] < ZERO)
{
wt1 = -ONE/Sdi[1];
usewt1 = TRUE;
}
else if (Sdi[1] > ZERO)
{
usewt1 = FALSE;
}
else
{
*ierr1 = 1100;
ermsg( "DSFIT", *ierr1, 0, "Require SD(1) .ne. Zero", '.' );
return;
}
/* Test ordering of XI() array.
* */
if (Xi[nxy] >= Xi[1])
{
direct = TRUE;
}
else
{
direct = FALSE;
nxyp1 = nxy + 1;
}
*ierr1 = 0;
for (i = 2; i <= nxy; i++)
{
if (direct)
{
if (Xi[i - 1] > Xi[i])
{
*ierr1 = 400;
goto L_50;
}
}
else
{
if (Xi[i - 1] < Xi[i])
{
*ierr1 = 400;
goto L_50;
}
}
}
L_50:
;
if (*ierr1 != 0)
{
ierm1( "DSFIT", *ierr1, 0, "Require abcissas, X(), to be sorted."
, "I", i, ',' );
ierv1( "NXY", nxy, ',' );
derv1( " X(1)", Xi[1], ',' );
derv1( "X(I-1)", Xi[i - 1], ',' );
derv1( " X(I)", Xi[i], ',' );
derv1( "X(NXY)", Xi[nxy], '.' );
goto L_200;
}
/* Begin loop to form equations for least-squares spline fit.
* */
irnow = 1;
k = 1;
ksize = 0;
iseg = 1;
left = iseg + korder - 1;
for (jpt = 1; jpt <= nxy; jpt++)
{
if (direct)
{
jpoint = jpt;
}
else
{
jpoint = nxyp1 - jpt;
}
if (k > ldw)
{
dbacc( w, ldw, nband, &irnow, ksize, iseg, &jtprev, &ierr2 );
if (ierr2 != 0)
{
*ierr1 = 700 + ierr2;
goto L_100;
}
if (irnow > ldw)
{
*ierr1 = 600;
ierm1( "DSFIT", *ierr1, 0, "Require LDW .ge. NC+2"
, "LDW", ldw, ',' );
ierv1( "NC", nc, '.' );
goto L_200;
}
k = irnow;
ksize = 0;
}
/* DO WHILE(XI(JPOINT) .ge. TKNOTS(LEFT+1) .and. ISEG .lt. ISGMAX) */
L_60:
if (Xi[jpoint] >= Tknots[left + 1] && iseg < isgmax)
{
dbacc( w, ldw, nband, &irnow, ksize, iseg, &jtprev, &ierr2 );
if (ierr2 != 0)
{
*ierr1 = 800 + ierr2;
goto L_100;
}
ksize = 0;
k = irnow;
iseg += 1;
left += 1;
goto L_60;
}
/* END WHILE
*
* Build one equation */
dsbasd( korder, left, tknots, Xi[jpoint], 0, p );
if (usewt1)
{
wt = wt1;
}
else
{
sdijp = Sdi[jpoint];
if (sdijp > ZERO)
{
wt = ONE/sdijp;
}
else
{
*ierr1 = 1200;
ermsg( "DSFIT", *ierr1, 0, "With SD(1) > 0 require all SD(I) > 0."
, ',' );
derv1( "SD(1)", Sdi[1], ',' );
ierv1( "I", jpoint, ',' );
derv1( "SD(I)", sdijp, '.' );
return;
}
}
for (j = 1; j <= korder; j++)
{
W(j - 1,k - 1) = P[j]*wt;
}
W(kordp1 - 1,k - 1) = Yi[jpoint]*wt;
/* End of build one equation */
k += 1;
ksize += 1;
}
dbacc( w, ldw, nband, &irnow, ksize, iseg, &jtprev, &ierr2 );
if (ierr2 != 0)
{
*ierr1 = 900 + ierr2;
goto L_100;
}
/* ALL DATA POINTS HAVE BEEN PROCESSED. CALL FOR SOLUTION.
* */
dbsol( 1, w, ldw, nband, irnow, jtprev, bcoef, nc, &rnorm, &ierr3 );
if (ierr3 != 0)
{
*ierr1 = 1000 + ierr2;
ermsg( "DSFIT", *ierr1, 0, "Singularity noted in DBSOL.",
'.' );
goto L_200;
}
dof = max( 1, nxy - nc );
*sigfac = rnorm/sqrt( dof );
return;
/* ERROR IN _BACC */
L_100:
ermsg( "DSFIT", *ierr1, 0, "Error detected in subroutine DBACC"
, '.' );
/* ERROR RETURN */
L_200:
for (i = 1; i <= nc; i++)
{
Bcoef[i] = ZERO;
}
return;
#undef W
} /* end of function */