/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:32:07 */
/*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 "dsfitc.h"
#include <float.h>
#include <stdlib.h>
#include <string.h>
/* 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*/ dsfitc(
byte ccode[][5],
double xi[],
double yi[],
double sd[],
long korder,
long ncoef,
double tknots[],
double bcoef[],
double *rnorm,
long iset[],
long info[],
double 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;
double 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 */
double *const Basis = &basis[0] - 1;
double *const Bcoef = &bcoef[0] - 1;
long *const Info = &info[0] - 1;
long *const Iset = &iset[0] - 1;
double *const Sd = &sd[0] - 1;
double *const Tknots = &tknots[0] - 1;
double *const Work = &work[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.
*>> 1996-03-30 DSFITC Krogh Added external statement.
*>> 1996-01-23 DSGITC Krogh Changes to simplify conversion to C.
*>> 1995-11-21 DSFITC Krogh Converted from SFTRAN to Fortran 77.
*>> 1994-11-16 CLL Add loops to zero debug arrays XT() and RT().
*>> 1994-10-19 DSFITC Krogh Changes to use M77CON
*>> 1994-01-31 DSFITC 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 DBLSE. 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 DBLSE. 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.
* ------------------------------------------------------------------
*--D 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( "DSFITC", ierr5, 0, "Require NCOEF .ge. 1", "NCOEF"
, ncoef, '.' );
goto L_500;
}
if (korder > KMAX)
{
ierr5 = 150;
ierm1( "DSFITC", 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( "DSFITC", ierr5, 0, "Require knots, TKNOTS(I), to be nondecreasing."
, "I", i, ',' );
derv1( "TKNOTS(I)", Tknots[i], ',' );
derv1( "TKNOTS(I+1)", Tknots[i + 1], '.' );
goto L_500;
}
}
for (i = 1; i <= ncoef; i++)
{
if (Tknots[i] >= Tknots[i + korder])
{
ierr5 = 250;
ierm1( "DSFITC", ierr5, 0, "Require TKNOTS(I) < TKNOTS(I+KORDER)."
, "I", i, ',' );
derv1( "TKNOTS(I)", Tknots[i], ',' );
derv1( "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( "DSFITC", 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( "DSFITC", 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( "DSFITC", 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( "DSFITC", 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( "DSFITC", ierr5, 0, "Quitting on bad values in CCODE()"
, '.' );
goto L_500;
}
goto L_60;
L_180:
;
/* end forever
* */
if (nbadcc != 0)
{
ermsg( "DSFITC", ierr5, 0, "Quitting on bad values in CCODE()"
, '.' );
goto L_500;
}
mtot = m1 + mfit;
ntot = ncoef + ns;
/* print*,'DSFITC.. 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( "DSFITC", 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 DSFIND.
* J1,J2 = Arbitrary starting values for use by DSBASI.
* */
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( "DSFITC", ierr5, 0, "With SD(1) > 0 require all SD(I) > 0."
, ',' );
derv1( "SD(1)", Sd[1], ',' );
ierv1( "I", ic, ',' );
derv1( "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 */
dsbasi( 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 */
dsfind( tknots, korder, ncoef + 1, x, &left, &mode );
dsbasd( korder, left, tknots, x, deriv, basis );
jcol = left - korder + 1;
/* print*,'DSFITC.. X =',X
* print*,'DSFITC.. 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)','DSFITC.. MTOT =',MTOT,', NTOT =',NTOT
* print'(1x/1x,a/1x)','DSFITC.. Matrix going to DBLSE:'
* 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 = pow(DBL_EPSILON,0.75e0);
dblse( work, mtot, mtot, ntot, m1, &Work[mn + 1], (double(*)[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( "DSFITC", ierr5, 0, "Error noted in DBLSE.", '.' );
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 */