/*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 "sc2fit.h"
#include <stdlib.h>
/* PARAMETER translations */
#define NBAND 4
#define NBAND1 (NBAND + 1)
#define ONE 1.0e0
#define ZERO 0.0e0
/* end of PARAMETER translations */
void /*FUNCTION*/ sc2fit(
float xi[],
float yi[],
float sdi[],
long nxy,
float b[],
long nb,
float *w,
long nw,
float yknot[],
float ypknot[],
float *sigfac,
long *ierr1)
{
#define W(I_,J_) (*(w+(I_)*(nw)+(J_)))
LOGICAL32 newseg, usewt1;
long int i, ierr2, ierr3, irnow, iseg, j, jpoint, jtprev, k, ksize,
n1, nparam;
float dof, p[4], rnorm, sdijp, wt, wt1;
/* OFFSET Vectors w/subscript range: 1 to dimension */
float *const B = &b[0] - 1;
float *const P = &p[0] - 1;
float *const Sdi = &sdi[0] - 1;
float *const Xi = &xi[0] - 1;
float *const Yi = &yi[0] - 1;
float *const Yknot = &yknot[0] - 1;
float *const Ypknot = &ypknot[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-01 SC2FIT Krogh Dim. SDI(*) instead NXY.
*>> 1995-11-21 SC2FIT Krogh Converted from SFTRAN to Fortran 77.
*>> 1994-10-19 SC2FIT Krogh Changes to use M77CON
*>> 1994-01-31 SC2FIT CLL Added test for SDI(i) .le. 0 when SDI(1) > 0.
*>> 1990-01-23 CLL Deleted ref to unused variable NX in call to IERM1
*>> 1989-10-20 CLL
*>> 1987-10-22 SC2FIT Lawson Initial code.
* Least squares fit to discrete data by a C-2 cubic spline.
* ------------------------------------------------------------------
* Algorithm and program designed by C.L.Lawson and R.J.Hanson.
* The general approach but not the complete code is given in
* 'SOLVING LEAST SQUARES PROBLEMS', by Lawson and Hanson,
* publ by Prentice-Hall, 1974.
* Programming and later changes and corrections by Lawson,Hanson,
* T.Lang, and D.Campbell, Sept 1968, Nov 1969, and Aug 1970.
* Modified 1968 Sept 17 to provide C-2 continuity.
* 1974 5/21, C.L.Lawson, Fixed bug that caused ISEG to get too big.
* Also changed to exit immidiately if B() 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.
* 1989-10-20 CLL Changed code so there are no RETURN statements in
* Sftran procedures. Previous code had such a RETURN that led to
* a warning diagnostic from a Cray compiler due to unreachable
* CONTINUE statements. Also introduced "c--" lines for CHGTYP.
* ------------------------------------------------------------------
* SUBROUTINE ARGUMENTS
*
* (XI(i),i=1,NXY) [in] Abcissas of data to be fitted. Require
* this data be ordered so X(i) .le. X(i+1).
*
* (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.
*
* (B(j),j=1,NB) [in] Breakpoints for the spline function,
* including endpoints. These breakpoints must be
* strictly increasing: B(j) .lt. B(j+1).
* It is required that all abcissas, XI(i), lie in the
* closed interval, [B(1), B(NB)].
*
* NB [in] No. of breakpoints, including endpoints.
* The no. of parameters in the least squares problem
* will be NB + 2.
* To have a nonsingular problem one must have
* NXY .ge. NB + 2, and the distribution of the breakpoints
* 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.
*
* W() [scratch] Work space dimensioned W(NW,5).
*
* NW [in] First dimension of W(). Must satisfy NW .ge. NB + 4
* Let KMAX denote the max no. of data abcissas, XI(i),
* in any one breakpoint interval, i.e. between B(j) and
* B(j+1) for some j. The subr will be more efficient
* if NW is at least NB + 3 + KMAX.
*
* (YKNOT(k) and YPKNOT(k),k=1,NB) [out] The subr will return
* values defining the fitted C2 spline curve in these arrays.
* These values and first derivatives of the fitted curve at the
* knot abcissae. YKNOT(j) = f(B(j)) and
* YPKNOT(j) = fprime(B(j)) for j = 1,...,NB.
* The user can then evaluate the fitted curve at any point by
* Hermite interpolation. See subrs DHINT or SHINT.
*
* 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 - (NB+2))
*
* IERR1 [out] Error status indicator. Note that IERR2 comes from
* SBACC and IERR3 comes from SBSOL.
*
* = 0 means no errors detected.
* = 100 means NB .lt. 2 .or. NXY .lt. NB+2
* = 200 means B(I) .ge. B(I+1)
* = 300 means NW .lt. NB+4
* = 400 means XI(I-1) .gt. XI(I)
* = 500 means B(1) .gt. XI(1) .or. B(NB) .lt. XI(NXY)
* = 600 means Need larger dimension NW.
* = 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 SBSOL. */
/* = 1100 means SDI(1) = zero.
* = 1200 means SDI(1) > zero and SDI(i) .le. zero for some i.
* ------------------------------------------------------------------
* Important internal variables.
*
* ISEG Index of current spline segment, starting with 1 for
* the first segment.
* 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.
* KSIZE Size of current block.
* JPOINT Current data pointer.
* ------------------------------------------------------------------
*--S replaces "?": ?C2FIT, ?BACC, ?C2BAS, ?ERM1, ?ERV1, ?BSOL, ?TRC2C
* Both versions use ERMSG, IERM1, IERV1
* Lower level subrs needed: (D/S)HTCC, (D/S)NRM2, ERFIN
* ------------------------------------------------------------------ */
/* ------------------------------------------------------------------ */
*ierr1 = 0;
/* EXIT IMMEDIATELY IF NB .lt. 2 OR NXY .LT NB+2 OR IF THE
* BREAKPOINTS ARE NOT STRICTLY INCREASING.
* */
nparam = nb + 2;
if (nb < 2 || nxy < nparam)
{
*ierr1 = 100;
ierm1( "SC2FIT", *ierr1, 0, "Require NB .ge. 2 and NXY .ge. NB+2"
, "NB", nb, ',' );
ierv1( "NXY", nxy, '.' );
goto L_300;
}
n1 = nb - 1;
for (i = 1; i <= n1; i++)
{
if (B[i] >= B[i + 1])
{
*ierr1 = 200;
ierm1( "SC2FIT", *ierr1, 0, "Require knots, B(I), to be strictly increasing."
, "I", i, ',' );
serv1( "B(I)", B[i], ',' );
serv1( "B(I+1)", B[i + 1], '.' );
goto L_300;
}
}
/* Require NW .ge. NB+4
* */
if (nw < nb + 4)
{
*ierr1 = 300;
ierm1( "SC2FIT", *ierr1, 0, "Require NW .ge. NB+4", "NW",
nw, ',' );
ierv1( "NB", nb, '.' );
goto L_300;
}
/* ------------------------------------------------------------------
* TEST SDI(1) */
if (Sdi[1] < ZERO)
{
wt1 = -ONE/Sdi[1];
usewt1 = TRUE;
}
else if (Sdi[1] > ZERO)
{
usewt1 = FALSE;
}
else
{
*ierr1 = 1100;
ermsg( "SC2FIT", *ierr1, 0, "Require SD(1) .ne. Zero", '.' );
return;
}
/* Test ordering of XI() array.
* */
for (i = 2; i <= nxy; i++)
{
if (Xi[i - 1] > Xi[i])
{
*ierr1 = 400;
ierm1( "SC2FIT", *ierr1, 0, "Require abcissas, X(I), to be nondecreasing."
, "I", i, ',' );
serv1( "X(I-1)", Xi[i - 1], ',' );
serv1( "X(I)", Xi[i], '.' );
goto L_300;
}
}
/* TEST THE FIRST AND LAST BREAKPOINT
* FOR BRACKETING THE DATA ABCISSAS.
* */
if (B[1] > Xi[1] || B[nb] < Xi[nxy])
{
*ierr1 = 500;
serm1( "SC2FIT", *ierr1, 0, "Require B(1) .LE. X(1) and B(NB) .ge. XI(NXY)"
, "B(1)", B[1], ',' );
serv1( "X(1)", Xi[1], ',' );
serv1( "B(NB)", B[nb], ',' );
serv1( "X(NXY)", Xi[nxy], '.' );
goto L_300;
}
/* BEGIN LOOP TO FORM EQUATIONS FOR C2 LEAST SQUARES FIT.
* */
irnow = 1;
k = 1;
ksize = 0;
newseg = TRUE;
iseg = 1;
for (jpoint = 1; jpoint <= nxy; jpoint++)
{
if (k > nw)
{
sbacc( w, nw, NBAND, &irnow, ksize, iseg, &jtprev, &ierr2 );
if (ierr2 != 0)
{
*ierr1 = 700 + ierr2;
goto L_200;
}
if (irnow > nw)
{
*ierr1 = 600;
ierm1( "SC2FIT", *ierr1, 0, "Need larger dimension NW."
, "NW", nw, '.' );
goto L_300;
}
k = irnow;
ksize = 0;
}
/* DO WHILE( XI(JPOINT) .gt. B(ISEG+1) ) */
L_80:
if (Xi[jpoint] > B[iseg + 1])
{
sbacc( w, nw, NBAND, &irnow, ksize, iseg, &jtprev, &ierr2 );
if (ierr2 != 0)
{
*ierr1 = 800 + ierr2;
goto L_200;
}
ksize = 0;
k = irnow;
iseg += 1;
newseg = TRUE;
goto L_80;
}
/* END WHILE
* Build one equation */
sc2bas( Xi[jpoint], 0, iseg, &newseg, b, nb, p );
if (usewt1)
{
wt = wt1;
}
else
{
sdijp = Sdi[jpoint];
if (sdijp > ZERO)
{
wt = ONE/sdijp;
}
else
{
*ierr1 = 1200;
ermsg( "SC2FIT", *ierr1, 0, "With SD(1) > 0 require all SD(I) > 0."
, ',' );
serv1( "SD(1)", Sdi[1], '.' );
ierv1( "I", jpoint, ',' );
serv1( "SD(I)", sdijp, '.' );
return;
}
}
for (j = 1; j <= 4; j++)
{
W(j - 1,k - 1) = P[j]*wt;
}
W(4,k - 1) = Yi[jpoint]*wt;
/* End of build one equation */
k += 1;
ksize += 1;
}
sbacc( w, nw, NBAND, &irnow, ksize, iseg, &jtprev, &ierr2 );
if (ierr2 != 0)
{
*ierr1 = 900 + ierr2;
goto L_200;
}
/* ALL DATA POINTS HAVE BEEN PROCESSED. CALL FOR SOLUTION.
* */
sbsol( 1, w, nw, NBAND, irnow, jtprev, &W(NBAND1 - 1,0), nparam,
&rnorm, &ierr3 );
if (ierr3 != 0)
{
*ierr1 = 1000 + ierr2;
ermsg( "SC2FIT", *ierr1, 0, "Singularity noted in SBSOL.",
'.' );
goto L_300;
}
dof = max( 1, nxy - nparam );
*sigfac = rnorm/sqrtf( dof );
/* TRANSFORM PARAMETERS TO Y,YPRIME BASIS
* */
strc2c( b, nb, &W(4,0), yknot, ypknot );
return;
/* Error in _BACC */
L_200:
ermsg( "SC2FIT", *ierr1, 0, "Error detected in subroutine SBACC"
, '.' );
/* Set YKNOT() & YPKNOT() to zero */
L_300:
for (i = 1; i <= nb; i++)
{
Yknot[i] = ZERO;
Ypknot[i] = ZERO;
}
return;
#undef W
} /* end of function */