/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:32:08 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "ssfitc.h" #include #include #include /* PARAMETER translations */ #define ATAB " AaNn!" #define DTAB "0123456789" #define ELIMIT 9 #define KINFO 7 #define KMAX 20 #define NVTAB "1234" #define ONE 1.0e0 #define RTAB "~=<>" #define UNBND 99.0e0 #define ZERO 0.0e0 /* end of PARAMETER translations */ void /*FUNCTION*/ ssfitc( byte ccode[][5], float xi[], float yi[], float sd[], long korder, long ncoef, float tknots[], float bcoef[], float *rnorm, long iset[], long info[], float work[]) { #define CCODE(I_,J_) (&ccode[I_][J_]) char _c0[2]; LOGICAL32 usewt1; long int _l0, active, deriv, fac, i, ic, ierr4, ierr5, ircon, irls, irow, iwbnd, iwbvec, iwcc, iwdiff, iwindx, iwjsta, iwrt, iwsiz, iwss, iwtnrm, iwwrk, iwxs, iwxt, iwz, j, j1, j2, jcol, js, kind, kprint, left, m1, mfit, minmn, mn, mode, mtot, nbadcc, need1, need2, ninfo, ns, nsetp, nt, ntot, nwork, relop; float basis[KMAX], rtval, sdic, tol, wt, wt1, x; static char msg1[20] = "CCODE(I)(1:1) = "; static char msg3[20] = "CCODE(I)(3:3) = "; static char msg4[20] = "CCODE(I)(4:4) = "; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Basis = &basis[0] - 1; float *const Bcoef = &bcoef[0] - 1; long *const Info = &info[0] - 1; long *const Iset = &iset[0] - 1; float *const Sd = &sd[0] - 1; float *const Tknots = &tknots[0] - 1; float *const Work = &work[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. *>> 1996-03-30 SSFITC Krogh Added external statement. *>> 1996-01-23 DSGITC Krogh Changes to simplify conversion to C. *>> 1995-11-21 SSFITC Krogh Converted from SFTRAN to Fortran 77. *>> 1994-11-16 CLL Add loops to zero debug arrays XT() and RT(). *>> 1994-10-19 SSFITC Krogh Changes to use M77CON *>> 1994-01-31 SSFITC CLL Added test for SD(i) .le. 0 when SD(1) > 0. *>> 1992-12-16 CLL Corrected formula for NEED2 and comments re WORK(). *>> 1992-11-12 C. L. Lawson, JPL Initializing LEFT, J1, J2. *>> 1992-10-27 C. L. Lawson, JPL *>> 1989-03-02 C. L. Lawson, JPL *>> 1989-02-23 C. L. Lawson, JPL *>> 1988-04-01 C. L. Lawson, JPL * * Weighted least squares fit to discrete data by a polynomial * spline function of order KORDER. The user can specify equality or * inequality constraints. The fitting equations as well as the * constraints may involve the value, a derivative of specified * order, or a definite integral of the spline function. A fitting * equation or constraint may involve function or derivative values * at different points, and relations between derivatives of * different orders. * * The order of a polynomial spline function is one * greater than the degree of its polynomial pieces. * Example: KORDER = 4 specifies a cubic spline function. * * The "proper fitting interval" is from A = TKNOTS(KORDER) to * B = TKNOTS(NT+1-KORDER). Extrapolation outside this interval * is permitted, but one must expect diminished accuracy at * extrapolated points. The given data or constraint * abcissas may be outside [A, B], and subsequent evaluations can be * done at points X outside [A, B]. * ------------------------------------------------------------------ * Specification of fitting and constraint equations. * * Let F denote the polynomial spline to be determined. Let the * quadruple (CCODE(i), X(i), Y(i), SD(i)) be called the ith * specification row. For each desired fitting equation or * constraint equation, (either of which we call a relation) the user * must specify one or two (consecutive) specification rows. * * CCODE(i) consists of 4 characters, that we call * KIND, DERIV, RELOP, and ACTIVE. * * KIND may be '1', '2', '3', or '4'. KIND determines the kind of * relation being specified. When KIND = '1' or '2' all * information for the relation is given in a single specification * row. We call this Row i in the following discussion. * When KIND = '3' or '4', two consecutive specification rows are use * We call these Rows i and i+1. We will complete the explanation of * KIND after describing DERIV, RELOP, and ACTIVE. * * DERIV may be '0', '1', ..., or '9'. This selects the order of * derivative of F to appear in the relation. '0' denotes the value * of F itself. * * RELOP may be '~', '=', '<', or '>'. * RELOP = '~' means the relation is to be a least-squares fitting * equation, '=' means an equality constraint equation, '<' means a * less-than-or-equal constraint equation, and '>' means a * greater-than-or-equal constraint equation. When RELOP = '~', * SD(i) specifies the a priori standard deviation of the right-side * member of the equation. When RELOP indicates a constraint * equation, SD(i) will be ignored. * * ACTIVE may be 'A', 'N', or '!'. Lower case 'a' and 'n' are * also accepted. ACTIVE = 'A' means the current specification row * is active, i.e., a relation will be generated from these * specifications. ACTIVE = 'N' means the specifications are not * active, i.e., no relation will be generated. ACTIVE = '!' means * the current row is inactive and there are no following rows, * i.e., this marks the end of the specification data. The user must * provide this termination signal. * When KIND = '3', or '4', meaning two specification rows are to * be interpreted together, the setting of ACTIVE = 'A' or * ACTIVE = 'N' must be consistent in these two rows. * Setting ACTIVE = '!' in only the first row is permitted, however, * since then the second row will not be accessed anyway. * * The forms of the relations selected by KIND are: * * KIND = 1: G(X(i)) RELOP Y(i) * * where G is the derivative of F selected by DERIV. * * KIND = 2: G(X(i)) - G(Y(i)) RELOP 0 * * where G is the derivative of F selected by DERIV. * * KIND = 3: G(X(i)) - Y(i+1) * H(X(i+1)) RELOP Y(i) * * where G is the derivative of F selected by DERIV(i) and * and H is the derivative of F selected by DERIV(i+1). In this * case KIND(i+1) and RELOP(i+1) are not used. * * KIND = 4: Integral from X(i) to X(i+1) of F RELOP Y(i) * */ /* In this case DERIV(i), KIND(i+1), DERIV(i+1), RELOP(i+1), and * Y(i+1) are not used. * ------------------------------------------------------------------ * The linear algebra methods were designed by C.L.Lawson and * R.J.Hanson. The method of representing spline functions is due * to Carl de Boor. References: * "SOLVING LEAST SQUARES PROBLEMS", by Lawson and Hanson, * Prentice-Hall, 1974. * "A PRACTICAL GUIDE TO SPLINES" by Carl de Boor, * Springer-Verlag, 1978. * The functionality and user interface of this subprogram are * modeled on the "French Curve" subroutine, FC, developed by Hanson * and Lawson at JPL in 1970, and the subsequent version developed by * Hanson at Sandia in 1979. * March 1988, CLL, JPL. Revised to conform to the Fortran 77 * standard. Intended for inclusion in the JPL MATH77 math library. * Feb 1989, CLL, JPL. Revised to use KIND = 3 to specify an * integral, and added new kind of relation using KIND = 4. * ------------------------------------------------------------------ * SUBROUTINE ARGUMENTS * * CCODE() [in, char*4] CCODE(i) is regarded as consisting of four * single-character fields. * CCODE(i)(1:1) = KIND = '1', '2', '3', '4'. * CCODE(i)(2:2) = DERIV = '0', '1', ..., '9'. * CCODE(i)(3:3) = RELOP = '~', '=', '<', '>'. * CCODE(i)(4:4) = ACTIVE = 'A', 'N', '!' * Where alphabetic characters are shown, the corresponding * lower case character is also acceptable. * * X() [in] Abcissas for specification of fitting or * constraint equations. * * Y() [in] Values or abcissas for specification * of fitting or constraint equations. * * SD() [in] Specifies the a priori standard deviation of error in the * right-side value in each fitting equation. * The weighted fitting algorithm will take account of these. * Optionally, the user may set SD(1) to a negative value. * Then this subr will use abs(SD(1)) as the standard deviation * for the right-side value in each fitting equation. In this * latter case the SD() array can be dimensioned SD(1). * Note that a negative value in SD(1) will always be interpreted * in this way by this subr, even if the associated RELOP is not * '~' or if ACTIVE is not 'A'. * * KORDER [in] Order of the spline basis functions. 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. * * NCOEF [in] No. of B-spline coefficients to be determined. * * (TKNOTS(j),j=1,NT, where NT = NCOEF+KORDER) [in] This is the deBoor * knot sequence for definition of the spline basis functions. * These values must be nondecreasing. * Repeated values are permitted, but values at * an index spacing of KORDER must be strictly increasing. * 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 and * convenient 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. * * BCOEF() [out] An array of length NCOEF into which the computed * coefficients defining the fitted curve will be stored. These * are coeffients relative to B-spline basis functions. * For I = 1, ..., NCOEF, the coefficient BCOEF(I) * is associated with the basis function whose support interval * runs from TKNOTS(I) to TKNOTS(I+KORDER). * * RNORM [out] RNORM := sqrt( sum over the indices i for which * fitting was requested of [( (yfit(i) - Y(i))/SD(i))**2]) * * ISET() [in integer] Array of length 3. * ISET(1) = NINFO, the dimension of INFO(). A sufficiently * large value for NINFO is 7 + 2*(NCOEF + NS). * ISET(2) = NWORK, the dimension of WORK(). A sufficiently * large value for NWORK can be computed as follows: * See definition of NS, M1, and MFIT below under INFO(). * NTOT = NCOEF + NS * MTOT = M1 + MFIT * MINMN = min(MTOT, NTOT) * NWORK = MTOT*NTOT + 3*MTOT + 6*NTOT + 3*MINMN + M1 * ISET(3) = KPRINT, a print flag in the range [0, 4]. It is * passed on to SBLSE. Larger values produce more printing. * * INFO() [out and scratch integer] The first 7 elements of INFO() * are used to return information about the problem. The * following 2*(NCOEF+NS) locations are used as scratch. * INFO(1) = IERR5, the Error status indicator. * Note that IERR4 comes from SBLSE. Possible values of IERR5 are * as follows: * */ /* = 0 means no errors detected. * = 100 means NCOEF .lt. 1 * = 200 means TKNOTS(I) .gt. TKNOTS(I+1) * = 250 means TKNOTS(I) .ge. TKNOTS(I+KORDER) * = 300 means NINFO or NWORK is too small. * = 500 means DERIV has bad value. * = 600 means RELOP has bad value. * = 700 means KIND has bad value. * = 800 means ACTIVE has bad value. * = 1000 + IERR4 means IERR4 .ne. 0 due to error * detected in _BLSE. * = 1100 means SD(1) = zero. * = 1200 means SD(1) > zero and SD(i) .le. zero for some i. * * INFO(2) = NEED1, the dimension needed for INFO(). * INFO(3) = NEED2, the dimension needed for WORK(). * INFO(4) = M1, the number of constraints rows in the matrix * representation of the problem. This will be a count of * the number of nonignored instances of CCODE(i) having * RELOP = '=', '<', or '>, and ACTIVE = 'A'. * INFO(5) = MFIT, the number of least-squares equations. * This will be a count of * the number of nonignored instances of CCODE(i) having * RELOP = '~' and ACTIVE = 'A'. * INFO(6) = NS, the number of slack variables. This will be a * count of the number of nonignored instances of CCODE(i) * having RELOP = '<' or '>, and ACTIVE = 'A'. * INFO(7) = NSETP, the number of variables in Set P at * termination. These variables are at values determined by * solution of a system of equations. The other NCOEF + NS * - NSETP variables will be at fixed values, either at one of * their bounds or at zero. * * WORK() [scratch] Work space dimensioned NWORK. * ------------------------------------------------------------------ * Important internal variables. * * BASIS() Array in which values of KORDER basis functions or their * will be stored. Dimensioned using the parameter KMAX. * This puts an upper limit on permissible KORDER. * JCOL Column of matrix into which first element of current * set of basis function values will be placed. * JCOL = LEFT - KORDER + 1. * KINFO Parameter specifying the number of locations at the * beginning of INFO() used for specific items of returned * information. Space beyond these locations is used for * scratch. * KMAX Intermal dimensioning parameter. The input value of KORDER * must not exceed KMAX. * KORDP1 = KORDER+1 * KSIZE Number of rows in current block. * LEFT Index of current spline segment. LEFT will satisfy * KORDER .le. LEFT .le. NCOEF. * The knot interval associated with index LEFT is from * T(LEFT) to T(LEFT+1). * Note that the union of these * segments is the "proper fitting interval". * ELIMIT Limit on number of errors in initial scan of specs before * quitting. * ------------------------------------------------------------------ *--S replaces "?": ?SFITC, ?BLSE, ?SBASI, ?SBASD, ?SFIND, ?ERV1 * Both versions use ERMOR, ERMSG, IERM1, IERV1 * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ ninfo = Iset[1]; nwork = Iset[2]; kprint = Iset[3]; nt = ncoef + korder; ierr5 = 0; /* Exit immediately if NCOEF .lt. 1 or if * the knots fail to be nondecreasing. * */ if (ncoef < 1) { ierr5 = 100; ierm1( "SSFITC", ierr5, 0, "Require NCOEF .ge. 1", "NCOEF" , ncoef, '.' ); goto L_500; } if (korder > KMAX) { ierr5 = 150; ierm1( "SSFITC", ierr5, 0, "Require KORDER .le. KMAX.", "KORDER" , korder, ',' ); ierv1( "KMAX", KMAX, '.' ); goto L_500; } for (i = 1; i <= (nt - 1); i++) { if (Tknots[i] > Tknots[i + 1]) { ierr5 = 200; ierm1( "SSFITC", ierr5, 0, "Require knots, TKNOTS(I), to be nondecreasing." , "I", i, ',' ); serv1( "TKNOTS(I)", Tknots[i], ',' ); serv1( "TKNOTS(I+1)", Tknots[i + 1], '.' ); goto L_500; } } for (i = 1; i <= ncoef; i++) { if (Tknots[i] >= Tknots[i + korder]) { ierr5 = 250; ierm1( "SSFITC", ierr5, 0, "Require TKNOTS(I) < TKNOTS(I+KORDER)." , "I", i, ',' ); serv1( "TKNOTS(I)", Tknots[i], ',' ); serv1( "TKNOTS(I+KORDER)", Tknots[i + korder], '.' ); goto L_500; } } /* ------------------------------------------------------------------ * TEST SD(1) */ if (Sd[1] < ZERO) { wt1 = -ONE/Sd[1]; usewt1 = TRUE; } else if (Sd[1] > ZERO) { usewt1 = FALSE; } else { ierr5 = 1100; ermsg( "SSFITC", ierr5, 0, "Require SD(1) .ne. Zero", '.' ); goto L_500; } /* . Determine M1, MFIT, NS, MTOT, and NTOT. * . M1 = number of non-ignored constraint specifications. * . MFIT = number of non-ignored least-squates equations. * . NS = number of non-ignored constraints that * . are inequality constraints and thus require a slack variable. * . MTOT = M1 + MFIT * . NTOT = NCOEF + NS * */ nbadcc = 0; m1 = 0; mfit = 0; ns = 0; i = 1; /* do forever */ L_60: ; kind = istrstr( NVTAB, STR1(_c0,CCODE(i - 1,0)[0]) ); deriv = istrstr( DTAB, STR1(_c0,CCODE(i - 1,0)[1]) ) - 1; relop = istrstr( RTAB, STR1(_c0,CCODE(i - 1,0)[2]) ); active = istrstr( ATAB, STR1(_c0,CCODE(i - 1,0)[3]) )/2; if (active == 3) goto L_180; /* . CCODE(I)(1:1) = 1 2 3 4 * . KIND = 1 2 3 4 * * . CCODE(I)(2:2) = 0 1 2 3 4 5 6 7 8 9 * . DERIV = 0 1 2 3 4 5 6 7 8 9 * * . CCODE(I)(3:3) = ~ = < > * . RELOP = 1 2 3 4 * * . CCODE(I)(4:4) = A a N n ! * . ACTIVE = 1 1 2 2 3 * */ if (active == 1) { /* do case(RELOP, 4) */ switch (relop) { case 1: goto L_110; case 2: goto L_120; case 3: goto L_130; case 4: goto L_140; } /* case other */ ierr5 = 600; ierm1( "SSFITC", ierr5, 0, "RELOP = CCODE(I)(3:3) has invalid value." , "I", i, ',' ); msg3[18] = CCODE(i - 1,0)[2]; ermor( msg3, '.' ); nbadcc += 1; goto L_150; /* case 1 */ L_110: mfit += 1; goto L_150; /* case 2 */ L_120: m1 += 1; goto L_150; /* case 3 */ L_130: m1 += 1; ns += 1; goto L_150; /* case 4 */ L_140: m1 += 1; ns += 1; L_150: ; /* end case */ } else if (active != 2) { ierr5 = 800; ierm1( "SSFITC", ierr5, 0, "ACTIVE = CCODE(I)(4:4) has invalid value." , "I", i, ',' ); msg4[18] = CCODE(i - 1,0)[3]; ermor( msg4, '.' ); nbadcc += 1; } if (kind == 1 || kind == 2) { i += 1; } else if (kind == 3 || kind == 4) { i += 2; } else { ierr5 = 700; ierm1( "SSFITC", ierr5, 0, "KIND = CCODE(I)(1:1) has invalid value." , "I", i, ',' ); msg1[18] = CCODE(i - 1,0)[0]; ermor( msg1, '.' ); nbadcc = ELIMIT + 1; } if (nbadcc > ELIMIT) { ermsg( "SSFITC", ierr5, 0, "Quitting on bad values in CCODE()" , '.' ); goto L_500; } goto L_60; L_180: ; /* end forever * */ if (nbadcc != 0) { ermsg( "SSFITC", ierr5, 0, "Quitting on bad values in CCODE()" , '.' ); goto L_500; } mtot = m1 + mfit; ntot = ncoef + ns; /* print*,'SSFITC.. M1, MFIT, NCOEF, NS =',M1, MFIT, NCOEF, NS */ minmn = min( mtot, ntot ); Info[4] = m1; Info[5] = mfit; Info[6] = ns; /* . Set indices to partition the work arrays INFO() and WORK(). * */ iwindx = 1 + KINFO; iwjsta = iwindx + ntot; need1 = iwjsta + ntot - 1; Info[2] = need1; mn = mtot*ntot; iwbvec = mn + 1; iwbnd = iwbvec + mtot; iwxs = iwbnd + 2*ntot; iwwrk = iwxs + ntot; iwsiz = iwwrk + ntot; iwtnrm = iwsiz + m1; iwz = iwtnrm + ntot; iwcc = iwz + minmn; iwss = iwcc + minmn; iwxt = iwss + minmn; iwrt = iwxt + ntot; iwdiff = iwrt + mtot; need2 = iwdiff + mtot - 1; Info[3] = need2; /* . Check NINFO and NWORK. * */ if (ninfo < need1 || nwork < need2) { ierr5 = 300; ierm1( "SSFITC", ierr5, 0, "Require NINFO .ge. NEED1 and NWORK .ge. NEED2." , "NINFO", ninfo, ',' ); ierv1( "NEED1", need1, ',' ); ierv1( "NWORK", nwork, ',' ); ierv1( "NEED2", need2, '.' ); goto L_500; } /* . Zero an MTOT x NTOT space for the matrix. */ for (i = 1; i <= mn; i++) { Work[i] = ZERO; } /* . Zero debug arrays XT() and RT() in WORK(). */ for (j = 1; j <= ntot; j++) { Work[iwxt + j - 1] = ZERO; } for (j = 1; j <= mtot; j++) { Work[iwrt + j - 1] = ZERO; } /* . Set bounds for variables. Storing into a 2 x NTOT space. * */ for (j = 1; j <= (2*ntot); j++) { Work[iwbnd - 1 + j] = UNBND; } for (j = 1; j <= ns; j++) { Work[iwbnd + (ncoef + j - 1)*2] = ZERO; } /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * Begin loop to form equations, both constraints and least-squares. * * Initialize JS = column index of previous slack variable. * IRCON = row index of previous constraint equation. * IRLS = row index of previous least-squares equation. * IC = index of current specification data. * LEFT = Arbitrary starting value for use by SSFIND. * J1,J2 = Arbitrary starting values for use by SSBASI. * */ js = ncoef; ircon = 0; irls = m1; ic = 1; j1 = 1; j2 = 1; left = 1; /* do forever */ L_300: ; active = istrstr( ATAB, STR1(_c0,CCODE(ic - 1,0)[3]) )/2; if (active == 3) goto L_400; if (active == 2) { ic += 1; goto L_300; } kind = istrstr( NVTAB, STR1(_c0,CCODE(ic - 1,0)[0]) ); relop = istrstr( RTAB, STR1(_c0,CCODE(ic - 1,0)[2]) ); /* . CCODE(I)(1:1) = 1 2 3 4 * . KIND = 1 2 3 4 * * . CCODE(I)(2:2) = 0 1 2 3 4 5 6 7 8 9 * . DERIV = 0 1 2 3 4 5 6 7 8 9 * * . CCODE(I)(3:3) = ~ = < > * . RELOP = 1 2 3 4 * * . CCODE(I)(4:4) = A a N n ! * . ACTIVE = 1 1 2 2 3 * * . Set matrix row index, IROW. * . Set weight, WT, if RELOP = 1. * . Store coefficient of +1 or -1 for slack variable if * . RELOP is 3 or 4. * */ if (relop == 1) { irls += 1; irow = irls; if (usewt1) { wt = wt1; } else { sdic = Sd[ic]; if (sdic > ZERO) { wt = ONE/sdic; } else { ierr5 = 1200; ermsg( "SSFITC", ierr5, 0, "With SD(1) > 0 require all SD(I) > 0." , ',' ); serv1( "SD(1)", Sd[1], ',' ); ierv1( "I", ic, ',' ); serv1( "SD(I)", sdic, '.' ); goto L_500; } } } else { ircon += 1; irow = ircon; if (relop == 3) { js += 1; Work[irow + (js - 1)*mtot] = ONE; } else if (relop == 4) { js += 1; Work[irow + (js - 1)*mtot] = -ONE; } } if (kind == 4) { /* Setup for integral */ ssbasi( korder, ncoef, tknots, Xi[ic], Xi[ic + 1], &j1, &j2, &Work[iwxs] ); for (j = j1; j <= j2; j++) { if (relop == 1) { Work[irow + (j - 1)*mtot] = Work[iwxs - 1 + j]*wt; } else { Work[irow + (j - 1)*mtot] = Work[iwxs - 1 + j]; } } if (relop == 1) { Work[mn + irow] = Yi[ic]*wt; } else { Work[mn + irow] = Yi[ic]; } /* End of setup for integral */ ic += 2; } else { /* Setup for value or derivative * . DERIV = CCODE()(2:2) = 0 1 2 3 4 5 6 7 8 9 * */ deriv = istrstr( DTAB, STR1(_c0,CCODE(ic - 1,0)[1]) ) - 1; x = Xi[ic]; fac = ONE; /* Negate KIND to flag that this is the first time accumulating values. */ kind = -kind; L_340: ; /* Accumulate values into matrix */ ssfind( tknots, korder, ncoef + 1, x, &left, &mode ); ssbasd( korder, left, tknots, x, deriv, basis ); jcol = left - korder + 1; /* print*,'SSFITC.. X =',X * print*,'SSFITC.. IC,LEFT,JCOL =',IC,LEFT,JCOL */ if (relop == 1) fac *= wt; for (j = 1; j <= korder; j++) { Work[irow + (jcol - 1)*mtot] += fac*Basis[j]; jcol += 1; } /* End of accumulate values into matrix */ if (kind < 0) { /* KIND is reset to positive here. */ kind = -kind; if (kind == 1) { rtval = Yi[ic]; } else { /* . Here KIND = 2 or 3 */ if (kind == 2) { fac = -ONE; rtval = ZERO; x = Yi[ic]; } else { deriv = istrstr( DTAB, STR1(_c0,CCODE(ic,0)[1]) ) - 1; fac = -Yi[ic + 1]; rtval = Yi[ic]; x = Xi[ic + 1]; } /* Go back to accumulate values a second time and last time on this iter. */ goto L_340; } } /* Set right side of relation. */ if (relop == 1) { Work[mn + irow] = rtval*wt; } else { Work[mn + irow] = rtval; } /* End of setup for value or derivative */ ic += 1; if (kind == 3) ic += 1; } goto L_300; L_400: ; /* end forever * . End loop to form equations. * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * print'(1x/1x,a,i5,a,i5)','SSFITC.. MTOT =',MTOT,', NTOT =',NTOT * print'(1x/1x,a/1x)','SSFITC.. Matrix going to SBLSE:' * do for I = 1,MTOT * print'(1x/1x,i5,3x,5g13.5/(9x,5g13.5))', * * I,(WORK(I+(J-1)*MTOT),J=1,NTOT+1) * end for * * . All points have been processed. Call for the solution. * */ tol = powf(FLT_EPSILON,0.75e0); sblse( work, mtot, mtot, ntot, m1, &Work[mn + 1], (float(*)[2])&Work[iwbnd], UNBND, kprint, tol, &ierr4, &Work[iwxs], rnorm, &nsetp, &Work[iwwrk], &Work[iwsiz], &Work[iwtnrm], &Work[iwz], &Work[iwcc], &Work[iwss], &Info[iwindx], &Info[iwjsta], &Work[iwxt], &Work[iwrt], &Work[iwdiff] ); /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ Info[7] = nsetp; if (ierr4 != 0) { ierr5 = 1000 + ierr4; ermsg( "SSFITC", ierr5, 0, "Error noted in SBLSE.", '.' ); goto L_500; } for (i = 1; i <= ncoef; i++) { Bcoef[i] = Work[iwxs - 1 + i]; } Info[1] = ierr5; return; /* ------------------------------------------------------------------ * ERROR RETURN */ L_500: for (i = 1; i <= ncoef; i++) { Bcoef[i] = ZERO; } Info[1] = ierr5; return; #undef CCODE } /* end of function */