/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:41 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "dilupm.h" #include #include /* PARAMETER translations */ #define LIDXSZ (MAXDIM*(2 + MAXDEG)) #define LTXTAB 9 #define LTXTAC 54 #define LTXTAD 100 #define LTXTAE 133 #define LTXTAF 180 #define LTXTAG 234 #define LTXTAH 289 #define LTXTAI 341 #define LTXTAJ 379 #define LTXTAK 440 #define LTXTAL 499 #define LTXTAM 550 #define MAXDEG 15 #define MAXDIM 10 #define MECONT 50 #define MEEMES 52 #define MEFVEC 61 #define MEIVEC 57 #define MEMDA1 27 #define MERET 51 #define METEXT 53 /* end of PARAMETER translations */ /* COMMON translations */ struct t_cdilup { double badpt, dx[MAXDEG + 1-(0)+1], ebnd, ebndr; long int kaos, kextrp, kgoc, lderiv, lexerr, lexit, ndegec, idx[MAXDEG + 1-(0)+1]; LOGICAL32 geterr; } cdilup; /* end of COMMON translations */ void /*FUNCTION*/ dilupm( long ndim, double x[], double *y, long ntab[], double xt[], double yt[], long ndeg[], long lup[], long iopt[], double eopt[]) { LOGICAL32 geterp; long int i, id, idxbas, ieeopt, iiflg, imaxd, iop2n, iopti[3], ip[MAXDIM], iprag, irag, ix, ixt, iybas, iydsav, iysav, iyt, j, k, kc[MAXDIM], kdim, kl, l, lckend, lertmp, li, licinf, lidx[LIDXSZ], linfo[MAXDIM-(0)+1], linfo0[MAXDIM-(0)+1], lkgoc[MAXDIM], lndeg[MAXDIM], lopt, lt, lupc, mtotd, mxd, ndimi, nt, ntabm, ntabxt, numd[MAXDIM], setext, _i, _r; static long int intchk[20-(0)+1]; double errdat[2], errsav[MAXDIM]; static char mtxtaa[3][202]={"DILUPM$BEstimated error = $F, Requested error = $F.$ENeed NDIM in interval [1, $I] but NDIM = $I.$ELUP($I) = $I; it should be < 4.$ERagged table badly specified in dimension $I.$EStart$ of ragged table", ", NTAB($I) = $I; it must be 0.$EMiddle of ragged table, NTAB($I) = $I; it must be >0.$EEnd of ragged table, NTAB($I) = $I; it must be -1.$EIOPT($I) = $I is not a valid option.$EIn dimension $I, first X", "T = XT($I) = last XT = XT($I) = $F.$EA total of only $I derivatives allowed, but $I requested.$ENeed derivative order in [0, $I], but it is = $I.$EPrevious error from DILUP$ occurred in dimension $I.$E"}; static char mtxtab[1][13]={"LUP(1:$M):$B"}; static char mtxtac[1][14]={"NDEG(1:$M):$B"}; static char mtxtad[1][14]={"NTAB(1:$M):$B"}; static char mtxtae[1][14]={"IOPT(1:$M):$B"}; static char mtxtaf[1][12]={"X((1:$M):$B"}; static long mloc[12]={LTXTAB,LTXTAC,LTXTAD,LTXTAE,LTXTAF,LTXTAG, LTXTAH,LTXTAI,LTXTAJ,LTXTAK,LTXTAL,LTXTAM}; static long mact[12]={MEEMES,0,0,0,MECONT,MEMDA1,0,METEXT,MEIVEC, 0,MECONT,MERET}; static long id4 = 4; static int _aini = 1; /* EQUIVALENCE translations */ long _e0[12]; long int *const kptd = (long*)_e0; long int *const stod = (long*)_e0; /* end of EQUIVALENCE translations */ /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Eopt = &eopt[0] - 1; double *const Errdat = &errdat[0] - 1; double *const Errsav = &errsav[0] - 1; long *const Iopt = &iopt[0] - 1; long *const Iopti = &iopti[0] - 1; long *const Ip = &ip[0] - 1; long *const Kc = &kc[0] - 1; long *const Lidx = &lidx[0] - 1; long *const Lkgoc = &lkgoc[0] - 1; long *const Lndeg = &lndeg[0] - 1; long *const Lup = &lup[0] - 1; long *const Mact = &mact[0] - 1; long *const Mloc = &mloc[0] - 1; long *const Ndeg = &ndeg[0] - 1; long *const Ntab = &ntab[0] - 1; long *const Numd = &numd[0] - 1; double *const X = &x[0] - 1; double *const Xt = &xt[0] - 1; double *const Yt = &yt[0] - 1; /* end of OFFSET VECTORS */ if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */ intchk[0] = 120; _aini = 0; } /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>> 2006-03-28 DILUPM Krogh NDEG odd was giving bad error est. *>> 1999-07-17 DILUPM Krogh Changed MAXDIM to 10. *>> 1997-09-30 DILUPM Krogh Fixed bug on setting NT for ragged case. *>> 1995-12-13 DILUPM Krogh Fixed bad Y value when getting derivative. *>> 1995-12-04 DILUPM Krogh Set so DILUP errors keep stop/print level. *>> 1995-11-10 DILUPM Krogh Fixed so char. data at col. 72 is not ' '. *>> 1994-11-11 DILUPM Krogh Declared all vars. *>> 1994-10-20 DILUPM Krogh Changes to use M77CON *>> 1994-09-12 DILUPM Krogh Added CHGTYP code. *>> 1993-04-28 DILUPM Krogh Additions for Conversion to C. *>> 1992-04-08 DILUPM Krogh Removed unused label 350, add ID4. *>> 1991-10-17 DILUPM Krogh Initial Code. * *--D replaces "?": ?ILUP, ?ILUPM, C?ILUP, ?MESS * * Multidimensional Polynomial Interpolation with Look Up * Design/Code by Fred T. Krogh, Jet Propulsion Laboratory, Pasadena, CA * Last Revision: March 1, 1991 * * This subroutine works by making multiple calls to DILUP. See comments * there for the kind of look up options and the kind of interpolation * options available. * * ***************** Formal Arguments *************************** * * If the data is not defined on a grid of values of the independent * variables, the table is said to be ragged. In order to describe the * situation for ragged tables, define * r = index of first NTAB(i) such that NTAB(i) < 0. * e = -NTAB(r) if r is defined, else e = 1000. * If e = 1000, data is on a grid and one can ignore most of what is said * in connection with ragged tables. In the case of ragged tables, * either r = NDIM, or for i .ge. r, NTAB(i) = -(i-1). * Below we refer to lexicographic order as the order in which certain * data is stored. Let i.1, i.2, ... i.NDIM be indices associated with * dimensions 1, to NDIM. Data, D, which depends on these indices, is * said to be in lexicographic order if for all valid values of the * indices the data is stored in consecutive location in memory, and data * for j.1, ..., j.NDIM precedes data for i.1, ..., i.NDIM if j.1 < i.1, * and in the case of equality if j.2 < i.2, etc. Thus for fixed values * of the first NDIM-1 indices, successive values of the last index will * be in successive memory locations. * * NDIM Number of dimensions for the independent variable (1 3*NDIM + 1 + space * needed to store information on the raggedness. Information for * each dimension with tabular information that depends on indices * selected for earlier dimensions starts with a 0, then the * number of table entries for the first value of the indices from * lower dimensions, then the number of entries for the next value * of the indices, etc. Here "first" and "next" means with the * indices ordered lexicographically, that is indices from the * lowest dimension vary first, and those from the highest * dimension last. Thus if (j,i) represent indices from the * second and the first dimension, the order is (1,1), (1,2), ..., * (2,1), (2,2), .... Information on raggedness is supplied first * for the first ragged dimension, then for the next, etc. * The 0 starting information for one dimension must follow * immediately after the last item from the previous dimension. * A -1 must follow the information on raggedness for the last * dimension. This is used to enhance the error checks. * XT An array giving values of the independent variable where Y is * is known. The organization of XT is defined by NTAB. * If the values of Y are known on a grid of points, this * organization is simple. * YT Array of dependent variable values. In the case of a grid, one * can think of YT as an NDIM dimensional array with YT(I(NDIM), * I(NDIM-1), ..., I(1)) being the value of Y at X = (XT(I(1)), * XT(I(2)), ..., XT(I(NDIM))). In the case of general ragged * tables (as for a grid too), the YT values are stored in * lexicographic order. * NDEG NDEG(i) defines the nominal degree of the polynomial to be used * in the ith dimension, with definition just like that for NDEG * in DILUP. * LUP LUP(i) defines the type of look up method for the ith dimension * exactly as LUP does in DILUP. * IOPT IOPT(1) is used to return a status. Values are as for DILUP, * with the addition of * - 1 Error estimate is > that requested error. * - 2 Bad value for NDIM. * - 3 Bad value for LUP(i). * - 4 Bad specification for a ragged table. * - 5 Ragged table does not start with a 0. * - 6 Bad value inside a ragged table. * - 7 Bad value at end of ragged table. */ /* - 8 Bad option index. * - 9 In some dimension, the first and last XT values are equal. * -10 Too many derivatives requested. * -11 Bad value for number of derivatives in some dimension. * -12 Problem with storage use in EOPT. * <-20 Had an error in DILUP. The value is the value returned * by DILUP - 20. * * IOPT(2) gives the dimension of EOPT. In addition to the first * location and the locations in EOPT used for options, DILUPM * needs a contiguous block of storage in EOPT of length > * 2*(2*NDIM - 1 + sum of NDEG(i), i = 1 to NDIM - 1). If * derivatives are being computed even more space is needed as is * described with option 3 below. * To get a message printed of the space required, set IOPT(2)=0. * * IOPT(k), k > 3 is used for the specifiation of options. * Options are as follows: * 0 End of the option list, must be the last item in IOPT. * 1 An error estimate is to be returned in EOPT(1) * 2 (Argument K2) K2 is a vector of length NDIM whose i-th entry * gives the polynomial degree to use when extrapolating with * respect to the i-th dimension. The default for the i-th entry * is 2 * max(1, NDEG(i)/2)). * 3 (Arguments K3, L3, M3) Let D(i.1, i.2, ..., i.NDIM) denote the * result of differentiating the interpolating polynomial i.j * times with respect to i.j, j = 1, 2, ..., NDIM. The i.j * satisfy the restriction: sum (i.j) .le. L3, and 0 .le. i.j .le. * M3.j, where M3 is a vector of length NDIM, 0 .le. M3.j .le. * min(NDEG(j), L3). Subject to these restrictions, the D's are * stored in lexicographic order of i.NDIM, ..., i.1 starting in * EOPT(K3). Thus for example if NDIM = 2, L3=2, and * M3.1 = M3.2 = 1, then EOPT(K3) = D(1, 0), EOPT(K3+1) = D(0, 1), * and EOPT(K3+2) = D(1, 1). If L3 had been 1, nothing would be * stored in EOPT(K3+2). * Extra storage in EOPT is required when this option is used, as * mentioned in the description of IOPT(2). The amount needed is * difficult to state for the general case. Let n.j denote the * total number of derivatives that must be computed in dimension * j in order to get the final derivatives in dimension 1. A * recurrence for computing the n.j is given in the subroutine * write-up. The extra storage required is = sum for j = 2, NDIM * of (NDEG(j-1) + 2) * n.j * If one does not want to deal with the complexities of figuring * out the storage, set K3 = 0, and a suggested organization will * be printed out and the program stopped. Or, set K3 = -1, and * the program will use the location it computes in the case K3 = * 0, and return the value used in IOPT where K3 was stored * initially. * 4 (Argument K4) The absolute and relative errors expected in YT * entries are specified in EOPT(K4) and EOPT(K4+1) respectively. * The values provided here are used in estimating the error in * the interpolation. An error estimate is returned in EOPT(1). * 5 (Argument K5) Do the interpolation to the accuracy requested * by the absolute error tolerance specified in EOPT(K5). An * attempt is made to keep the final error < EOPT(K5). Standard * polynomial interpolation is done, but here NDEG() gives the * maximum polynomial degree to use in the interpolation. If * EOPT(K5) .le. 0., IOPT(1) is not set to -1, and no error * message is generated due to an unsatisfied accuracy request. * An error estimate is returned in EOPT(1). * 6 (Argument K6) Do not use a point in the interpolation if the * corresponding value of YT = EOPT(K6). * EOPT Array used to return an error estimate, and used for options. * EOPT(1) if an error estimate is returned, this contains a crude * estimate of the error in the interpolation. * * ***************** Variables Referenced *********************** * * ABS Fortran intrinsic -- absolute value. * BADPT In common CDILUP, see documentation in DILUP. * DX In common CDILUP, see documentation for DILUP. * EBND In common CDILUP, see documentation for DILUP. * EBNDR In common CDILUP, see documentation for DILUP. * EOPT Formal argument, see above. * ERRDAT Place to store floating point for error messages. * ERRSAV Place for saving error estimates. * GETERP Usual value for GETERR. * GETERR In common CDILUP, see documentation for DILUP. * I Temporary index. * ID Current place to pick up derivatives from higher dimension in * order to interpolate derivative in a lower dimension. * ID4 = 4 defined in data. To avoid passing literal to DILUP. * IDX In common CDILUP, see documentation for DILUP. * IDXBAS Base location for storing values of DX in EOPT. ( =IOPT(2)+1) * The first DX is stored in EOPT(IDXBAS+1). * IEEOPT Index for absolute and relative error in EOPT. * IIFLG Value of flag to be returned in IOPT(1). * IMAXD IOPT(IMAXD+KDIM) gives the maximum derivative to be computed * with respect to dimension KDIM. * INTCHK Array used for checking use of optional space in EOPT. * IOP2N Value to store in IOPTI(2) when interpolating in the last * dimension. = 255 if bad points possible, else = 254. * IOPT Formal argument, see above. * IOPTI Internal array used to pass options to DILUP. * IP Array used to compute current indices into XT and YT. * IPRAG Part of a term used in computing IP when KDIM>IRAG, and * NTAB(KDIM) > 0. * IRAG = index of first ragged table, = 1000 if table not ragged. * IX Current point index w.r.t. the current dimension. */ /* IXT Pointer to XT information for the current dimension. * IYBAS One less than first location for storing values interpolated * for y in EOPT. * IYDSAV As for IYSAV below, but for saving a derivative value. * IYSAV Location in EOPT for saving interpolated value. * IYT Location in EOPT where vector of interpolated Y's are stored. * J Index of a 0 in NTAB which is followed by user information * specifying the number of data points as a function of * previously used values from XT. * K Temporary index. * KAOS In common CDILUP. Keeps track of state on variable order. * = 0 No variable order -- this value only set in DILUPM. * = 1 First time * = 2 Second time * = 3 Just had an increase or a very weak decrease in the error. * = 4 Just had a strong decrease in the error. * On an error with IIFLG .lt. 0, in DILUP this is set to -10 * stop * level - the print level. * KC Array used for counters for computing derivatives. * KDIM Index of dimension currently active. * KEXTRP In common CDILUP, see documentation in DILUP. * KGOC In common CDILUP, see documentation for DILUP. * KL Largest total derivative as starting from the highest dimension * and working to dimension 1 in getting storage needed for * derivatives. * KPTD KPTD(KDIM) = pointer to extra derivative information for * dimension KDIM. KPTD(KDIM-1) is place to store deriv. * values for dim. KDIM. * L Temporary index. * LCKEND End of space used for checking storage in EOPT. * LDERIV In common CDILUP, see documentation for DILUP. * LERTMP Location in EOPT where temp. error info. is stored (2 values). * LEXERR In common CDILUP, see documentation for DILUP. * LEXIT In common CDILUP. Defines action after finding location in XT. * 0 Don't compute anything for Y, just return. (Used for DILUPM.) * 1 Usual case, just interpolate. * >1 Compute LEXIT-1 derivatives of the interpolant. * LI Last location examined in IOPT(). * LICINF Base for indexing information that depends on the current * point index in the current dimension. * LIDX Array containing indices of the points selected for all * dimensions up to the current one. Also used in computing * KPTD. * LIDXSZ Parameter giving the amount of space allocated for LIDX. * LINFO LINFO(KDIM) contains a pointer to the index currently being * worked on in dimension KDIM. * LINFO0 LINFO0(KDIM) gives value to use for LICINF when in dim. KDIM, = * 1 + sum of (2 + NDEG(J)) for J = 1, 2, ..., KDIM-1. * LKGOC LKGOC(KDIM) contains value of KGOC for interp. in dim. KDIM. * LNDEG LNDEG(KDIM) contains value of NDEG for interp. in dim. KDIM. * LOPT Value from current location in IOPT. * LT As much of NTAB as is known to be safe for printing error info. * LTXTxx Parameter names of this form were generated by PMESS in making * up text for error messages and are used to locate various parts * of the message. * LUP Formal argument, see above. * LUPC Value passed to DILUP for LUP, the look up method. * MACT Array giving actions for printing error messages, see MESS. * MAX Fortran intrinsic -- maximum * MAXDEG Parameter -- Gives maximum degree of polynomial interpolation. * MAXDIM Paramter giving the maximum number of dimensions allowed. * MECONT Parameter telling MESS an error message is to be continued. * MEEMES Parameter telling MESS to print an error message. * MEIVEC Parameter telling MESS to print an integer vector. * MEFVEC Parameter telling MESS to print an floating point vector. * MEMDA1 Parameter telling MESS that next item is integer data to be * made available for output. * MERET Parameter telling MESS that this ends the error message. * MESS Subroutine for printing error messages. * METEXT Parameter telling MESS that data from MTXTAA is to printed. * MLOC Array giving locations of start of text for error messages. * MTOTD Maximum order of total derivative to be computed. * MTXTAx Character arrays holding error message text for MESS. * MXD Upper bound on number of derivatives to compute when computing * derivatives based on information from higher dimensions. * NDEG Formal argument, see above. * NDEGEC In common CDILUP, see documentation for DILUP. * NDIM Formal argument, see above. * NDIMI Internal value for NDIM = number of dimensions. * NT Number of points for the current call to DILUP. * NTAB Formal argument, see above. * NTABM = NTABXT + NDIMI -- NTAB(NTABM+I) is 1 + the index of the last * word in NTAB required for ragged table storage through dim. I. * NTABXT = NDIMI+1 -- NTAB(NTABXT) = index of first ragged table, =1000 * if there is no ragged table. NTABN(NTABXT+I) = base address * accessing XT information. * NUMD NUMD(j) gives the number of different derivatives for * dimension j. * OPTCHK Subroutine for checking on space allocation for options. * SETEXT If not 0, gives IOPT(SETEXT+KDIM) gives the value to use for * KEXTRP in dimension KDIM. * DILUP One dimensional interpolation subroutine. * DMESS Calls MESS and prints floating data in error messages. * STOD Index in EOPT where derivatives are to be stored. * X Formal argument, see above. * XT Formal argument, see above. * Y Formal argument, see above. * YT Formal argument, see above. * * ************* Formal Variable Declarations *************** */ /* ************* Common Block and Parameter ***************** * * MAXDEG must be odd and > 2. */ /* ************* Local Variables **************************** * */ /* ************************ Error Message Stuff and Data **************** * * Parameter defined below are all defined in the error message program * MESS and DMESS. * */ /* ********* Error message text *************** *[Last 2 letters of Param. name] [Text generating message.] *AA DILUPM$B *AB Estimated error = $F, Requested error = $F.$E *AC Need NDIM in interval [1, $I] but NDIM = $I.$E *AD LUP($I) = $I; it should be < 4.$E *AE Ragged table badly specified in dimension $I.$E *AF Start of ragged table, NTAB($I) = $I; it must be 0.$E *AG Middle of ragged table, NTAB($I) = $I; it must be >0.$E *AH End of ragged table, NTAB($I) = $I; it must be -1.$E *AI IOPT($I) = $I is not a valid option.$E *AJ In dimension $I, first XT = XT($I) = last XT = XT($I) = $F.$E *AK A total of only $I derivatives allowed, but $I requested.$E *AL Need derivative order in [0, $I], but it is = $I.$E *AM Previous error from DILUP occurred in dimension $I.$E * $ *AN LUP(1:$M):$B * $ *AO NDEG(1:$M):$B * $ *AP NTAB(1:$M):$B * $ *AQ IOPT(1:$M):$B * $ *AR X((1:$M):$B */ /* 1 2 3 4 5 6 7 8 9 10 */ /* ***************** Start of executable code ******************* * */ li = 2; ndimi = ndim; ntabxt = ndimi + 1; ntabm = ntabxt + ndimi; lt = ntabm; if (Ntab[ntabxt] <= 0) { /* Fix NTAB on the first call. */ if (ndimi > MAXDIM) { intchk[1] = MAXDIM; intchk[2] = ndimi; iiflg = -2; goto L_400; } Ntab[ntabxt] = 1000; Ntab[ntabxt + 1] = 1; for (i = 1; i <= ndimi; i++) { if (Lup[i] >= 4) { /* LUP(I) is out of range -- fatal error. */ intchk[1] = i; intchk[2] = Lup[i]; iiflg = -3; goto L_400; } l = Ntab[i]; if (l > 0) { /* Table is not ragged up to this point. */ if (i == ndimi) goto L_50; /* If table is ragged save pointer to ragged info. */ if (Ntab[ndimi] < 0) Ntab[ntabm + i] = ntabm + ndimi; /* Get pointer to start of XT data for next dimension */ if (Lup[i] == 3) { Ntab[ntabxt + i + 1] = Ntab[ntabxt + i] + 2; } else { Ntab[ntabxt + i + 1] = Ntab[ntabxt + i] + l; } } else if ((l == 0) || ((l != 1 - i) && ((i != ndimi) || (l <= -i)))) { /* Problem in specifying raggedness. */ lt = ndimi; intchk[1] = i; iiflg = -4; goto L_400; } else { j = Ntab[ntabm + i - 1]; if (Ntab[i - 1] > 0) { Ntab[ntabxt] = -Ntab[i]; /* Get K = number of NTAB entries of extra data */ k = Ntab[1]; for (l = 2; l <= (-l); l++) { k *= Ntab[l]; } } else { k = Ntab[j - 1]; } if (Ntab[j] != 0) { lt = j; intchk[1] = j; intchk[2] = Ntab[j]; iiflg = -5; goto L_400; } /* Change data to be the partial sum of the original data. */ lt = j + k; for (l = j + 1; l <= lt; l++) { if (Ntab[l] <= 0) { /* Error, bad value inside ragged table info. */ intchk[1] = l; intchk[2] = Ntab[l]; iiflg = -6; goto L_400; } Ntab[l] += Ntab[l - 1]; } if (i == ndimi) { lt += 1; if (Ntab[lt] == -1) goto L_50; /* Error -- No end tag where needed. */ intchk[1] = lt; intchk[2] = Ntab[lt]; iiflg = -7; goto L_400; } /* Save index of 0th NTAB entry for extra data for next dim. */ Ntab[ntabm + i] = j + k + 1; /* Get pointer to start of XT data for next dimension */ if (Lup[i] == 3) { Ntab[ntabxt + i + 1] = Ntab[ntabxt + i] + 2*k; } else { Ntab[ntabxt + i + 1] = Ntab[ntabxt + i] + Ntab[j + k]; } } } } L_50: ; iiflg = 0; linfo[0] = 0; linfo0[0] = 0; linfo0[1] = 1; for (i = 2; i <= ndimi; i++) { linfo0[i] = linfo0[i - 1] + Ndeg[i - 1] + 2; } /* Initialize option flags */ geterp = FALSE; setext = 0; cdilup.kaos = 0; ieeopt = 0; *stod = 0; cdilup.lderiv = 0; iop2n = 254; irag = Ntab[ntabxt]; intchk[4] = 0; intchk[6] = 0; intchk[5] = -2*(linfo0[ndimi] + 1); lckend = 7; goto L_70; L_62: intchk[lckend + 2] = 1; L_64: intchk[lckend + 1] = Iopt[li]; L_65: intchk[lckend] = lopt; lckend += 3; L_70: ; li += 1; lopt = Iopt[li]; if (lopt != 0) { switch (lopt) { case 1: goto L_310; case 2: goto L_320; case 3: goto L_330; case 4: goto L_300; case 5: goto L_340; case 6: goto L_360; } intchk[1] = li; intchk[2] = lopt; iiflg = -8; goto L_400; } intchk[1] = lckend; intchk[2] = Iopt[2]; intchk[3] = 1; intchk[lckend] = 1; optchk( intchk, iopt, "DILUPM / EOPT$E" ); if (intchk[1] < 0) { Iopt[1] = -12; return; } j = lckend; /* Store values for storage locations determined by OPTCHK. */ L_75: j += 1; k = intchk[j]; if (intchk[k] == 0) { lertmp = intchk[intchk[lckend + 1] + 1]; } else { *stod = intchk[intchk[lckend + 2] + 1]; } if (j != intchk[lckend]) goto L_75; idxbas = lertmp + 1; iybas = idxbas + linfo0[ndimi]; if (*stod != 0) { k = iybas + linfo0[ndimi]; for (j = 1; j <= ndimi; j++) { kptd[j] += k; } cdilup.lderiv = kptd[ndimi - 1]; } iysav = iybas; for (i = lertmp; i <= iybas; i++) { /* Ensure in DILUP, YT will be defined. (It's ref. but not used.) */ Eopt[i] = 0; } Errsav[1] = 0.e0; iyt = idxbas; cdilup.geterr = geterp; kdim = 1; ixt = 1; lupc = Lup[1]; Ip[1] = 0; L_100: nt = Ntab[kdim]; L_110: cdilup.lexit = 0; /* Setup for variable order. */ if (cdilup.kaos != 0) { cdilup.kaos = 1; Iopti[3] = Ndeg[kdim]; } /* Set KEXTRP to indicate how to extrapolate */ cdilup.kextrp = -1; if (setext != 0) cdilup.kextrp = Iopt[setext + kdim]; if (kdim != ndimi) goto L_180; /* Interpolate in the last dimension. * Flag that want results on the exit */ cdilup.lexit = 1; /* Setup for getting derivatives. */ if (cdilup.lderiv > 0) cdilup.lexit = Iopt[imaxd + kdim] + 1; /* Indicate where error info. on Y is found (if any) */ cdilup.lexerr = ieeopt; /* Indicate whether there are bad points. */ Iopti[2] = iop2n; dilup( X[kdim], &Eopt[iysav], nt, &Xt[ixt], &Yt[Ip[kdim] + 1], Ndeg[kdim], &lupc, iopti, eopt ); if (lupc < 0) Lup[kdim] = lupc; if (cdilup.kaos < 0) goto L_390; iiflg = max( iiflg, Iopti[1] ); L_130: ; /* Back up to lower dimension */ kdim -= 1; if (kdim == 0) { *y = Eopt[iysav]; Iopt[1] = iiflg; if (Iopti[1] >= 0) return; iiflg = Iopti[1]; Errdat[1] = Eopt[1]; Errdat[2] = cdilup.ebnd; goto L_400; } if (geterp) Errsav[kdim] = fmax( Errsav[kdim], Eopt[1] ); cdilup.ndegec = Lndeg[kdim]; if (cdilup.kaos != 0) { if (linfo[kdim] >= linfo0[kdim] + 3) { cdilup.ndegec = linfo[kdim] - linfo0[kdim]; goto L_140; } } if (linfo[kdim] < linfo0[kdim] + cdilup.ndegec) goto L_220; /* Restore information in the common block */ L_140: cdilup.kgoc = labs( Lkgoc[kdim] ); licinf = linfo0[kdim]; linfo[kdim + 1] = linfo0[kdim + 1]; for (i = 0; i <= cdilup.ndegec; i++) { cdilup.idx[i] = Lidx[i + licinf]; cdilup.dx[i] = Eopt[i + idxbas + licinf]; } iysav = iybas + linfo[kdim - 1]; iyt = iybas + linfo0[kdim]; /* Flag that already have the Y's in order. */ lupc = 4; /* Flag that want results on the exit */ cdilup.lexit = 1; /* Setup for getting derivatives. */ if (cdilup.lderiv != 0) { if (cdilup.lderiv >= -kdim) { cdilup.lderiv = kptd[kdim - 1]; if (kdim != 1) cdilup.lderiv += Numd[kdim]*(linfo[kdim - 1] - linfo0[kdim - 1]); cdilup.lexit = Iopt[imaxd + kdim] + 1; } } /* Indicate where error info. on Y is found (if any) */ if (geterp) { cdilup.lexerr = lertmp; Eopt[cdilup.lexerr] = Errsav[kdim]; Eopt[cdilup.lexerr + 1] = 0.e0; } /* Setup for variable order. */ if (cdilup.kaos != 0) { cdilup.kaos = 1; Iopti[3] = cdilup.ndegec; } /* Indicate no bad points. */ L_180: Iopti[2] = 254; if ((Ndeg[kdim]%2) == 1) cdilup.ndegec = Ndeg[kdim]; dilup( X[kdim], &Eopt[iysav], nt, &Xt[ixt], &Eopt[iyt], Ndeg[kdim], &lupc, iopti, eopt ); if (cdilup.kaos < 0) goto L_390; iiflg = max( iiflg, Iopti[1] ); if (cdilup.lexit == 0) { /* Just selected points for this dimension */ if (lupc < 0) Lup[kdim] = lupc; /* Save info. from the common block. */ Lkgoc[kdim] = cdilup.kgoc; Lndeg[kdim] = cdilup.ndegec; licinf = linfo0[kdim]; linfo[kdim] = licinf; for (i = 0; i <= cdilup.ndegec; i++) { Lidx[i + licinf] = cdilup.idx[i]; Eopt[i + idxbas + licinf] = cdilup.dx[i]; } goto L_230; } else { /* Just got result for this dimension * IOPTI(2) = 0, if variable order and need more points. */ if (Iopti[2] == 0) goto L_220; if (cdilup.lderiv > 0) { cdilup.geterr = FALSE; id = kptd[kdim]; for (i = kdim + 1; i <= ndimi; i++) { Kc[i] = 0; } mxd = mtotd; L_214: iydsav = cdilup.lderiv + cdilup.lexit - 1; l = kdim + 1; L_215: if ((Kc[l] == Iopt[imaxd + l]) || (mxd == 0)) { mxd += Kc[l]; Kc[l] = 0; l += 1; if (l <= ndimi) goto L_215; goto L_218; } Kc[l] += 1; mxd -= 1; cdilup.lexit = min( mxd, Iopt[imaxd + kdim] ) + 1; cdilup.lderiv = iydsav + 1; for (i = 0; i <= cdilup.ndegec; i++) { Eopt[iyt + i] = Eopt[id + Numd[kdim + 1]*i]; } id += 1; cdilup.kgoc = labs( cdilup.kgoc ); if ((Ndeg[kdim]%2) == 1) cdilup.ndegec = Ndeg[kdim]; /* Interpolate in an outer dimension. */ dilup( X[kdim], &Eopt[iydsav], nt, &Xt[ixt], &Eopt[iyt], Ndeg[kdim], &id4, iopti, eopt ); goto L_214; L_218: cdilup.geterr = geterp; } goto L_130; } L_220: linfo[kdim] += 1; if (cdilup.lderiv > 0) { if (linfo[kdim] == linfo0[kdim] + cdilup.ndegec) { if (Lkgoc[kdim] == -2) { cdilup.lderiv = -kdim; goto L_230; } } if (kdim == ndimi - 1) { cdilup.lderiv += Numd[ndimi]; } else { cdilup.lderiv = kptd[ndimi - 1]; } } L_230: ix = Lidx[linfo[kdim]]; /* Get here after getting next point in this dimension * Compute IP() data for getting XT and YT indices */ if (kdim < irag) { Ip[kdim + 1] = Ntab[kdim + 1]*(ix - 1 + Ip[kdim]); } else if (kdim == irag) { k = Ntab[ntabm + kdim] + Ip[kdim] + ix; Ip[kdim + 1] = Ntab[k - 1]; iprag = Ntab[k] - Ntab[k - 1]; } else if (Ntab[kdim] > 0) { Ip[kdim + 1] = Ntab[kdim]*Ip[kdim] + (ix - 1)*iprag; } else { Ip[kdim + 1] = Ntab[Ntab[ntabm + kdim] + Ip[kdim] + ix - 1]; } iysav = iybas + linfo[kdim]; kdim += 1; Errsav[kdim] = 0.e0; if ((cdilup.lexit == 0) || (Ntab[kdim] < 0)) { ixt = Ntab[ntabxt + kdim]; lupc = Lup[kdim]; if (Ntab[kdim] >= 0) goto L_100; k = -Ntab[kdim]; k = Ntab[ntabm + k] + Ip[k]; nt = Ntab[k + ix] - Ntab[k + ix - 1]; if (lupc == 3) { ixt += 2*(Ip[-Ntab[kdim]] + Lidx[linfo[-Ntab[kdim]]] - 1); } else { ixt += Ip[1 - Ntab[kdim]]; } goto L_110; } if (kdim == ndim) { lupc = Lup[kdim]; goto L_110; } linfo[kdim] = linfo0[kdim]; goto L_230; /* Process various options * Specify error in Y. */ L_300: li += 1; ieeopt = Iopt[li]; intchk[lckend + 2] = 2; geterp = TRUE; goto L_64; /* Return an error estimate */ L_310: geterp = TRUE; goto L_70; /* Setup to set extrapolation degree */ L_320: setext = li; li += ndimi; goto L_70; /* Setup for computing extra derivatives */ L_330: *stod = Iopt[li + 1]; intchk[lckend + 1] = *stod; imaxd = li + 2; mtotd = Iopt[imaxd]; if (mtotd >= LIDXSZ) { intchk[1] = LIDXSZ; intchk[2] = mtotd; iiflg = -10; goto L_400; } kl = Iopt[imaxd + ndimi]; k = min( mtotd, Ndeg[ndimi] ); if ((kl < 0) || (kl > k)) { intchk[1] = k; intchk[2] = kl; iiflg = -11; goto L_400; } for (k = 0; k <= mtotd; k++) { Lidx[k + 1] = min( k, kl ) + 1; } Numd[ndimi] = kl; for (j = ndimi - 1; j >= 1; j--) { l = Iopt[imaxd + j]; k = min( mtotd, Ndeg[j] ); if ((l < 0) || (l > k)) { intchk[1] = k; intchk[2] = l; iiflg = -11; goto L_400; } for (k = mtotd; k >= (l + 1); k--) { Lidx[k + 1] -= Lidx[k - l]; } for (k = 1; k <= mtotd; k++) { Lidx[k + 1] += Lidx[k]; } kl = min( kl + l, mtotd ); Numd[j] = Lidx[kl + 1] - 1; } intchk[lckend + 2] = Numd[1]; if (*stod <= 0) { intchk[lckend + 1] = -Numd[1]; intchk[lckend + 2] = li + 1; } kptd[1] = 0; for (j = 2; j <= ndimi; j++) { kptd[j] = kptd[j - 1] + Numd[j]*(Ndeg[j - 1] + 2); } intchk[5] -= kptd[ndimi]; li = imaxd + ndimi; goto L_65; /* Setup for variable order. */ L_340: li += 1; cdilup.ebnd = Eopt[Iopt[li]]; cdilup.ebndr = 0.e0; cdilup.kaos = 1; geterp = TRUE; goto L_62; /* Setup to indicate bad points */ L_360: li += 1; cdilup.badpt = Eopt[Iopt[li]]; iop2n = 255; goto L_62; /* Error Processing */ L_390: iiflg = Iopti[1] - 20; intchk[1] = kdim; Mact[4] = Mloc[12]; Mact[2] = -cdilup.kaos; goto L_410; L_400: Mact[4] = Mloc[-iiflg]; Mact[2] = 88; if (iiflg == -1) Mact[2] = 25; L_410: Iopt[1] = iiflg; Mact[3] = -iiflg; dmess( mact, (char*)mtxtaa,202, &intchk[1], errdat ); if (iiflg < -2) { Mact[7] = ndimi; Mact[10] = ndimi; mess( &Mact[6], (char*)mtxtab,13, lup ); mess( &Mact[6], (char*)mtxtac,14, ndeg ); Mact[9] = MEFVEC; dmess( &Mact[6], (char*)mtxtaf,12, ndeg, x ); Mact[9] = MEIVEC; Mact[7] = lt; Mact[10] = lt; mess( &Mact[6], (char*)mtxtad,14, ntab ); if (li > 3) { Mact[7] = li; Mact[10] = li; mess( &Mact[6], (char*)mtxtae,14, iopt ); } } mess( &Mact[12], (char*)mtxtaa,202, &intchk[1] ); return; } /* end of function */