/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:59 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "ssfit.h" #include /* PARAMETER translations */ #define KMAX 20 #define ONE 1.0e0 #define ZERO 0.0e0 /* end of PARAMETER translations */ void /*FUNCTION*/ ssfit( float xi[], float yi[], float sdi[], long nxy, long korder, long nc, float tknots[], float bcoef[], float *sigfac, long *ierr1, long ldw, float *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; float dof, p[KMAX], rnorm, sdijp, wt, wt1; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Bcoef = &bcoef[0] - 1; float *const P = &p[0] - 1; float *const Sdi = &sdi[0] - 1; float *const Tknots = &tknots[0] - 1; float *const Xi = &xi[0] - 1; float *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 SSFIT Krogh Declared SDI(*) so ref. with SD(1) o.k. *>> 1995-11-21 SSFIT Krogh Converted from SFTRAN to Fortran 77. *>> 1994-10-19 SSFIT Krogh Changes to use M77CON *>> 1994-01-31 SSFIT CLL Add test for SDI(i) .le. 0 when SDI(1) > 0. *>> 1993-01-12 CLL Using SSBASD 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 * SBACC and IERR3 comes from SBSOL. * * = 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. * ------------------------------------------------------------------ *--S 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( "SSFIT", *ierr1, 0, "Require NC .ge. 1 and NXY .ge. NC" , "NC", nc, ',' ); ierv1( "NXY", nxy, '.' ); goto L_200; } if (korder > KMAX) { *ierr1 = 150; ierm1( "SSFIT", *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( "SSFIT", *ierr1, 0, "Require knots, TKNOTS(I), to be nondecreasing." , "I", i, ',' ); serv1( "TKNOTS(I)", Tknots[i], ',' ); serv1( "TKNOTS(I+1)", Tknots[i + 1], '.' ); goto L_200; } } for (i = 1; i <= nc; i++) { if (Tknots[i] >= Tknots[i + korder]) { *ierr1 = 250; ierm1( "SSFIT", *ierr1, 0, "Require TKNOTS(I) < TKNOTS(I+KORDER)." , "I", i, ',' ); serv1( "TKNOTS(I)", Tknots[i], ',' ); serv1( "TKNOTS(I+KORDER)", Tknots[i + korder], '.' ); goto L_200; } } /* Require LDW .ge. NC+2 * */ if (ldw < nc + 2) { *ierr1 = 300; ierm1( "SSFIT", *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( "SSFIT", *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( "SSFIT", *ierr1, 0, "Require abcissas, X(), to be sorted." , "I", i, ',' ); ierv1( "NXY", nxy, ',' ); serv1( " X(1)", Xi[1], ',' ); serv1( "X(I-1)", Xi[i - 1], ',' ); serv1( " X(I)", Xi[i], ',' ); serv1( "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) { sbacc( w, ldw, nband, &irnow, ksize, iseg, &jtprev, &ierr2 ); if (ierr2 != 0) { *ierr1 = 700 + ierr2; goto L_100; } if (irnow > ldw) { *ierr1 = 600; ierm1( "SSFIT", *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) { sbacc( 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 */ ssbasd( 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( "SSFIT", *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 <= 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; } sbacc( 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. * */ sbsol( 1, w, ldw, nband, irnow, jtprev, bcoef, nc, &rnorm, &ierr3 ); if (ierr3 != 0) { *ierr1 = 1000 + ierr2; ermsg( "SSFIT", *ierr1, 0, "Singularity noted in SBSOL.", '.' ); goto L_200; } dof = max( 1, nxy - nc ); *sigfac = rnorm/sqrtf( dof ); return; /* ERROR IN _BACC */ L_100: ermsg( "SSFIT", *ierr1, 0, "Error detected in subroutine SBACC" , '.' ); /* ERROR RETURN */ L_200: for (i = 1; i <= nc; i++) { Bcoef[i] = ZERO; } return; #undef W } /* end of function */