1 subroutine SILUPM (NDIM,X,Y,NTAB,XT,YT,NDEG,LUP,IOPT,EOPT) 2c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. 3c ALL RIGHTS RESERVED. 4c Based on Government Sponsored Research NAS7-03001. 5c>> 2006-03-28 SILUPM Krogh NDEG odd was giving bad error est. 6c>> 1999-07-17 SILUPM Krogh Changed MAXDIM to 10. 7c>> 1997-09-30 SILUPM Krogh Fixed bug on setting NT for ragged case. 8C>> 1995-12-13 SILUPM Krogh Fixed bad Y value when getting derivative. 9C>> 1995-12-04 SILUPM Krogh Set so SILUP errors keep stop/print level. 10C>> 1995-11-10 SILUPM Krogh Fixed so char. data at col. 72 is not ' '. 11C>> 1994-11-11 SILUPM Krogh Declared all vars. 12c>> 1994-10-20 SILUPM Krogh Changes to use M77CON 13c>> 1994-09-12 SILUPM Krogh Added CHGTYP code. 14c>> 1993-04-28 SILUPM Krogh Additions for Conversion to C. 15c>> 1992-04-08 SILUPM Krogh Removed unused label 350, add ID4. 16c>> 1991-10-17 SILUPM Krogh Initial Code. 17c 18c--S replaces "?": ?ILUP, ?ILUPM, C?ILUP, ?MESS 19c 20c Multidimensional Polynomial Interpolation with Look Up 21c Design/Code by Fred T. Krogh, Jet Propulsion Laboratory, Pasadena, CA 22c Last Revision: March 1, 1991 23c 24c This subroutine works by making multiple calls to SILUP. See comments 25c there for the kind of look up options and the kind of interpolation 26c options available. 27c 28c ***************** Formal Arguments *************************** 29c 30c If the data is not defined on a grid of values of the independent 31c variables, the table is said to be ragged. In order to describe the 32c situation for ragged tables, define 33c r = index of first NTAB(i) such that NTAB(i) < 0. 34c e = -NTAB(r) if r is defined, else e = 1000. 35c If e = 1000, data is on a grid and one can ignore most of what is said 36c in connection with ragged tables. In the case of ragged tables, 37c either r = NDIM, or for i .ge. r, NTAB(i) = -(i-1). 38c Below we refer to lexicographic order as the order in which certain 39c data is stored. Let i.1, i.2, ... i.NDIM be indices associated with 40c dimensions 1, to NDIM. Data, D, which depends on these indices, is 41c said to be in lexicographic order if for all valid values of the 42c indices the data is stored in consecutive location in memory, and data 43c for j.1, ..., j.NDIM precedes data for i.1, ..., i.NDIM if j.1 < i.1, 44c and in the case of equality if j.2 < i.2, etc. Thus for fixed values 45c of the first NDIM-1 indices, successive values of the last index will 46c be in successive memory locations. 47c 48c NDIM Number of dimensions for the independent variable (1<NDIM<10). 49c The upper limit is artificial to catch bad values. 50c X Vector of dimension at least NDIM giving the place where the 51c value of the interpolant is desired. 52c Y Value of the interpolant. The interpolating function is always 53c a piecewise polynomial of NDIM variables. 54c NTAB Defines organization of the table. In the case of a grid, 55c NTAB(1 to NDIM) contain n.1, ..., n.NDIM, the number of points 56c with respect to the various dimensions. NTAB(NDIM+1) must be 57c 0 initially. On the first call it is set to e as defined 58c above. Also on the first call, values of NTAB(NDIM+2:2*NDIM+1) 59c are set to values that are used to find the XT values to pass 60c to SILUP. In the case of ragged tables more space is required. 61c In this case NTAB must have dimension > 3*NDIM + 1 + space 62c needed to store information on the raggedness. Information for 63c each dimension with tabular information that depends on indices 64c selected for earlier dimensions starts with a 0, then the 65c number of table entries for the first value of the indices from 66c lower dimensions, then the number of entries for the next value 67c of the indices, etc. Here "first" and "next" means with the 68c indices ordered lexicographically, that is indices from the 69c lowest dimension vary first, and those from the highest 70c dimension last. Thus if (j,i) represent indices from the 71c second and the first dimension, the order is (1,1), (1,2), ..., 72c (2,1), (2,2), .... Information on raggedness is supplied first 73c for the first ragged dimension, then for the next, etc. 74c The 0 starting information for one dimension must follow 75c immediately after the last item from the previous dimension. 76c A -1 must follow the information on raggedness for the last 77c dimension. This is used to enhance the error checks. 78c XT An array giving values of the independent variable where Y is 79c is known. The organization of XT is defined by NTAB. 80c If the values of Y are known on a grid of points, this 81c organization is simple. 82c YT Array of dependent variable values. In the case of a grid, one 83c can think of YT as an NDIM dimensional array with YT(I(NDIM), 84c I(NDIM-1), ..., I(1)) being the value of Y at X = (XT(I(1)), 85c XT(I(2)), ..., XT(I(NDIM))). In the case of general ragged 86c tables (as for a grid too), the YT values are stored in 87c lexicographic order. 88c NDEG NDEG(i) defines the nominal degree of the polynomial to be used 89c in the ith dimension, with definition just like that for NDEG 90c in SILUP. 91c LUP LUP(i) defines the type of look up method for the ith dimension 92c exactly as LUP does in SILUP. 93c IOPT IOPT(1) is used to return a status. Values are as for SILUP, 94c with the addition of 95c - 1 Error estimate is > that requested error. 96c - 2 Bad value for NDIM. 97c - 3 Bad value for LUP(i). 98c - 4 Bad specification for a ragged table. 99c - 5 Ragged table does not start with a 0. 100c - 6 Bad value inside a ragged table. 101c - 7 Bad value at end of ragged table. 102c - 8 Bad option index. 103c - 9 In some dimension, the first and last XT values are equal. 104c -10 Too many derivatives requested. 105c -11 Bad value for number of derivatives in some dimension. 106c -12 Problem with storage use in EOPT. 107c <-20 Had an error in SILUP. The value is the value returned 108c by SILUP - 20. 109c 110c IOPT(2) gives the dimension of EOPT. In addition to the first 111c location and the locations in EOPT used for options, SILUPM 112c needs a contiguous block of storage in EOPT of length > 113c 2*(2*NDIM - 1 + sum of NDEG(i), i = 1 to NDIM - 1). If 114c derivatives are being computed even more space is needed as is 115c described with option 3 below. 116c To get a message printed of the space required, set IOPT(2)=0. 117c 118c IOPT(k), k > 3 is used for the specifiation of options. 119c Options are as follows: 120c 0 End of the option list, must be the last item in IOPT. 121c 1 An error estimate is to be returned in EOPT(1) 122c 2 (Argument K2) K2 is a vector of length NDIM whose i-th entry 123c gives the polynomial degree to use when extrapolating with 124c respect to the i-th dimension. The default for the i-th entry 125c is 2 * max(1, NDEG(i)/2)). 126c 3 (Arguments K3, L3, M3) Let D(i.1, i.2, ..., i.NDIM) denote the 127c result of differentiating the interpolating polynomial i.j 128c times with respect to i.j, j = 1, 2, ..., NDIM. The i.j 129c satisfy the restriction: sum (i.j) .le. L3, and 0 .le. i.j .le. 130c M3.j, where M3 is a vector of length NDIM, 0 .le. M3.j .le. 131c min(NDEG(j), L3). Subject to these restrictions, the D's are 132c stored in lexicographic order of i.NDIM, ..., i.1 starting in 133c EOPT(K3). Thus for example if NDIM = 2, L3=2, and 134c M3.1 = M3.2 = 1, then EOPT(K3) = D(1, 0), EOPT(K3+1) = D(0, 1), 135c and EOPT(K3+2) = D(1, 1). If L3 had been 1, nothing would be 136c stored in EOPT(K3+2). 137c Extra storage in EOPT is required when this option is used, as 138c mentioned in the description of IOPT(2). The amount needed is 139c difficult to state for the general case. Let n.j denote the 140c total number of derivatives that must be computed in dimension 141c j in order to get the final derivatives in dimension 1. A 142c recurrence for computing the n.j is given in the subroutine 143c write-up. The extra storage required is = sum for j = 2, NDIM 144c of (NDEG(j-1) + 2) * n.j 145c If one does not want to deal with the complexities of figuring 146c out the storage, set K3 = 0, and a suggested organization will 147c be printed out and the program stopped. Or, set K3 = -1, and 148c the program will use the location it computes in the case K3 = 149c 0, and return the value used in IOPT where K3 was stored 150c initially. 151c 4 (Argument K4) The absolute and relative errors expected in YT 152c entries are specified in EOPT(K4) and EOPT(K4+1) respectively. 153c The values provided here are used in estimating the error in 154c the interpolation. An error estimate is returned in EOPT(1). 155c 5 (Argument K5) Do the interpolation to the accuracy requested 156c by the absolute error tolerance specified in EOPT(K5). An 157c attempt is made to keep the final error < EOPT(K5). Standard 158c polynomial interpolation is done, but here NDEG() gives the 159c maximum polynomial degree to use in the interpolation. If 160c EOPT(K5) .le. 0., IOPT(1) is not set to -1, and no error 161c message is generated due to an unsatisfied accuracy request. 162c An error estimate is returned in EOPT(1). 163c 6 (Argument K6) Do not use a point in the interpolation if the 164c corresponding value of YT = EOPT(K6). 165c EOPT Array used to return an error estimate, and used for options. 166c EOPT(1) if an error estimate is returned, this contains a crude 167c estimate of the error in the interpolation. 168c 169c ***************** Variables Referenced *********************** 170c 171c ABS Fortran intrinsic -- absolute value. 172c BADPT In common CSILUP, see documentation in SILUP. 173c DX In common CSILUP, see documentation for SILUP. 174c EBND In common CSILUP, see documentation for SILUP. 175c EBNDR In common CSILUP, see documentation for SILUP. 176c EOPT Formal argument, see above. 177c ERRDAT Place to store floating point for error messages. 178c ERRSAV Place for saving error estimates. 179c GETERP Usual value for GETERR. 180c GETERR In common CSILUP, see documentation for SILUP. 181c I Temporary index. 182c ID Current place to pick up derivatives from higher dimension in 183c order to interpolate derivative in a lower dimension. 184c ID4 = 4 defined in data. To avoid passing literal to SILUP. 185c IDX In common CSILUP, see documentation for SILUP. 186c IDXBAS Base location for storing values of DX in EOPT. ( =IOPT(2)+1) 187c The first DX is stored in EOPT(IDXBAS+1). 188c IEEOPT Index for absolute and relative error in EOPT. 189c IIFLG Value of flag to be returned in IOPT(1). 190c IMAXD IOPT(IMAXD+KDIM) gives the maximum derivative to be computed 191c with respect to dimension KDIM. 192c INTCHK Array used for checking use of optional space in EOPT. 193c IOP2N Value to store in IOPTI(2) when interpolating in the last 194c dimension. = 255 if bad points possible, else = 254. 195c IOPT Formal argument, see above. 196c IOPTI Internal array used to pass options to SILUP. 197c IP Array used to compute current indices into XT and YT. 198c IPRAG Part of a term used in computing IP when KDIM>IRAG, and 199c NTAB(KDIM) > 0. 200c IRAG = index of first ragged table, = 1000 if table not ragged. 201c IX Current point index w.r.t. the current dimension. 202c IXT Pointer to XT information for the current dimension. 203c IYBAS One less than first location for storing values interpolated 204c for y in EOPT. 205c IYDSAV As for IYSAV below, but for saving a derivative value. 206c IYSAV Location in EOPT for saving interpolated value. 207c IYT Location in EOPT where vector of interpolated Y's are stored. 208c J Index of a 0 in NTAB which is followed by user information 209c specifying the number of data points as a function of 210c previously used values from XT. 211c K Temporary index. 212c KAOS In common CSILUP. Keeps track of state on variable order. 213c = 0 No variable order -- this value only set in SILUPM. 214c = 1 First time 215c = 2 Second time 216c = 3 Just had an increase or a very weak decrease in the error. 217c = 4 Just had a strong decrease in the error. 218c On an error with IIFLG .lt. 0, in SILUP this is set to -10 * stop 219c level - the print level. 220c KC Array used for counters for computing derivatives. 221c KDIM Index of dimension currently active. 222c KEXTRP In common CSILUP, see documentation in SILUP. 223c KGOC In common CSILUP, see documentation for SILUP. 224c KL Largest total derivative as starting from the highest dimension 225c and working to dimension 1 in getting storage needed for 226c derivatives. 227c KPTD KPTD(KDIM) = pointer to extra derivative information for 228c dimension KDIM. KPTD(KDIM-1) is place to store deriv. 229c values for dim. KDIM. 230c L Temporary index. 231c LCKEND End of space used for checking storage in EOPT. 232c LDERIV In common CSILUP, see documentation for SILUP. 233c LERTMP Location in EOPT where temp. error info. is stored (2 values). 234c LEXERR In common CSILUP, see documentation for SILUP. 235c LEXIT In common CSILUP. Defines action after finding location in XT. 236c 0 Don't compute anything for Y, just return. (Used for SILUPM.) 237c 1 Usual case, just interpolate. 238c >1 Compute LEXIT-1 derivatives of the interpolant. 239c LI Last location examined in IOPT(). 240c LICINF Base for indexing information that depends on the current 241c point index in the current dimension. 242c LIDX Array containing indices of the points selected for all 243c dimensions up to the current one. Also used in computing 244c KPTD. 245c LIDXSZ Parameter giving the amount of space allocated for LIDX. 246c LINFO LINFO(KDIM) contains a pointer to the index currently being 247c worked on in dimension KDIM. 248c LINFO0 LINFO0(KDIM) gives value to use for LICINF when in dim. KDIM, = 249c 1 + sum of (2 + NDEG(J)) for J = 1, 2, ..., KDIM-1. 250c LKGOC LKGOC(KDIM) contains value of KGOC for interp. in dim. KDIM. 251c LNDEG LNDEG(KDIM) contains value of NDEG for interp. in dim. KDIM. 252c LOPT Value from current location in IOPT. 253c LT As much of NTAB as is known to be safe for printing error info. 254c LTXTxx Parameter names of this form were generated by PMESS in making 255c up text for error messages and are used to locate various parts 256c of the message. 257c LUP Formal argument, see above. 258c LUPC Value passed to SILUP for LUP, the look up method. 259c MACT Array giving actions for printing error messages, see MESS. 260c MAX Fortran intrinsic -- maximum 261c MAXDEG Parameter -- Gives maximum degree of polynomial interpolation. 262c MAXDIM Paramter giving the maximum number of dimensions allowed. 263c MECONT Parameter telling MESS an error message is to be continued. 264c MEEMES Parameter telling MESS to print an error message. 265c MEIVEC Parameter telling MESS to print an integer vector. 266c MEFVEC Parameter telling MESS to print an floating point vector. 267c MEMDA1 Parameter telling MESS that next item is integer data to be 268c made available for output. 269c MERET Parameter telling MESS that this ends the error message. 270c MESS Subroutine for printing error messages. 271c METEXT Parameter telling MESS that data from MTXTAA is to printed. 272c MLOC Array giving locations of start of text for error messages. 273c MTOTD Maximum order of total derivative to be computed. 274c MTXTAx Character arrays holding error message text for MESS. 275c MXD Upper bound on number of derivatives to compute when computing 276c derivatives based on information from higher dimensions. 277c NDEG Formal argument, see above. 278c NDEGEC In common CSILUP, see documentation for SILUP. 279c NDIM Formal argument, see above. 280c NDIMI Internal value for NDIM = number of dimensions. 281c NT Number of points for the current call to SILUP. 282c NTAB Formal argument, see above. 283c NTABM = NTABXT + NDIMI -- NTAB(NTABM+I) is 1 + the index of the last 284c word in NTAB required for ragged table storage through dim. I. 285c NTABXT = NDIMI+1 -- NTAB(NTABXT) = index of first ragged table, =1000 286c if there is no ragged table. NTABN(NTABXT+I) = base address 287c accessing XT information. 288c NUMD NUMD(j) gives the number of different derivatives for 289c dimension j. 290c OPTCHK Subroutine for checking on space allocation for options. 291c SETEXT If not 0, gives IOPT(SETEXT+KDIM) gives the value to use for 292c KEXTRP in dimension KDIM. 293c SILUP One dimensional interpolation subroutine. 294c SMESS Calls MESS and prints floating data in error messages. 295c STOD Index in EOPT where derivatives are to be stored. 296c X Formal argument, see above. 297c XT Formal argument, see above. 298c Y Formal argument, see above. 299c YT Formal argument, see above. 300c 301c ************* Formal Variable Declarations *************** 302c 303 integer NDIM, NTAB(*), NDEG(*), LUP(*), IOPT(*) 304 real X(NDIM), Y, XT(*), YT(*), EOPT(*) 305c 306c ************* Common Block and Parameter ***************** 307c 308c MAXDEG must be odd and > 2. 309 integer MAXDEG 310 parameter (MAXDEG = 15) 311c 312 logical GETERR, GETERP 313 integer KAOS, KEXTRP, KGOC,LDERIV,LEXERR,LEXIT,NDEGEC, 314 1 IDX(0:MAXDEG+1) 315 real BADPT, DX(0:MAXDEG+1), EBND, EBNDR 316 common / CSILUP / BADPT, DX, EBND, EBNDR, KAOS, KEXTRP, KGOC, 317 1 LDERIV, LEXERR,LEXIT,NDEGEC,IDX,GETERR 318c 319c ************* Local Variables **************************** 320c 321 INTEGER MAXDIM, LIDXSZ 322 parameter (MAXDIM = 10) 323 parameter (LIDXSZ = MAXDIM*(2 + MAXDEG)) 324 integer IP(MAXDIM), LINFO(0:MAXDIM), 325 1 LINFO0(0:MAXDIM), LKGOC(MAXDIM), LNDEG(MAXDIM) 326 integer LIDX(LIDXSZ) 327 integer I, IDXBAS, IEEOPT, IIFLG, IMAXD, INTCHK(0:20), 328 1 IOP2N, IOPTI(3), IPRAG, IRAG, IX, IYBAS, IYDSAV, IYSAV, IYT, J, 329 2 K, KC(MAXDIM), KDIM, KPTD(0:MAXDIM), L, LERTMP, LI, LT, 330 3 LICINF, LOPT, LUPC, MTOTD, NDIMI, NT, NTABM, NTABXT, 331 4 NUMD(MAXDIM), SETEXT, STOD 332 integer ID4, LCKEND, IXT, ID, MXD, KL 333 equivalence (KPTD(0), STOD) 334 real ERRSAV(MAXDIM) 335c 336c ************************ Error Message Stuff and Data **************** 337c 338c Parameter defined below are all defined in the error message program 339c MESS and SMESS. 340c 341 integer MEMDA1, MECONT, MERET, MEEMES, METEXT, MEIVEC, MEFVEC 342 parameter (MEMDA1 =27) 343 parameter (MECONT =50) 344 parameter (MERET =51) 345 parameter (MEEMES =52) 346 parameter (METEXT =53) 347 parameter (MEIVEC =57) 348 parameter (MEFVEC =61) 349c 350 integer MLOC(12) 351 integer MACT(12) 352 real ERRDAT(2) 353c 354c ********* Error message text *************** 355c[Last 2 letters of Param. name] [Text generating message.] 356cAA SILUPM$B 357cAB Estimated error = $F, Requested error = $F.$E 358cAC Need NDIM in interval [1, $I] but NDIM = $I.$E 359cAD LUP($I) = $I; it should be < 4.$E 360cAE Ragged table badly specified in dimension $I.$E 361cAF Start of ragged table, NTAB($I) = $I; it must be 0.$E 362cAG Middle of ragged table, NTAB($I) = $I; it must be >0.$E 363cAH End of ragged table, NTAB($I) = $I; it must be -1.$E 364cAI IOPT($I) = $I is not a valid option.$E 365cAJ In dimension $I, first XT = XT($I) = last XT = XT($I) = $F.$E 366cAK A total of only $I derivatives allowed, but $I requested.$E 367cAL Need derivative order in [0, $I], but it is = $I.$E 368cAM Previous error from SILUP occurred in dimension $I.$E 369c $ 370cAN LUP(1:$M):$B 371c $ 372cAO NDEG(1:$M):$B 373c $ 374cAP NTAB(1:$M):$B 375c $ 376cAQ IOPT(1:$M):$B 377c $ 378cAR X((1:$M):$B 379 integer LTXTAA,LTXTAB,LTXTAC,LTXTAD,LTXTAE,LTXTAF,LTXTAG,LTXTAH, 380 * LTXTAI,LTXTAJ,LTXTAK,LTXTAL,LTXTAM,LTXTAN,LTXTAO,LTXTAP,LTXTAQ, 381 * LTXTAR 382 parameter (LTXTAA= 1,LTXTAB= 9,LTXTAC= 54,LTXTAD=100,LTXTAE=133, 383 * LTXTAF=180,LTXTAG=234,LTXTAH=289,LTXTAI=341,LTXTAJ=379, 384 * LTXTAK=440,LTXTAL=499,LTXTAM=550,LTXTAN= 1,LTXTAO= 1, 385 * LTXTAP= 1,LTXTAQ= 1,LTXTAR= 1) 386 character MTXTAA(3) * (201) 387 character MTXTAB(1) * (12) 388 character MTXTAC(1) * (13) 389 character MTXTAD(1) * (13) 390 character MTXTAE(1) * (13) 391 character MTXTAF(1) * (11) 392 data MTXTAA/'SILUPM$BEstimated error = $F, Requested error = $F.$E 393 *Need NDIM in interval [1, $I] but NDIM = $I.$ELUP($I) = $I; it sho 394 *uld be < 4.$ERagged table badly specified in dimension $I.$EStart$ 395 * of ragged table',', NTAB($I) = $I; it must be 0.$EMiddle of ragge 396 *d table, NTAB($I) = $I; it must be >0.$EEnd of ragged table, NTAB( 397 *$I) = $I; it must be -1.$EIOPT($I) = $I is not a valid option.$EIn 398 * dimension $I, first X','T = XT($I) = last XT = XT($I) = $F.$EA to 399 *tal of only $I derivatives allowed, but $I requested.$ENeed deriva 400 *tive order in [0, $I], but it is = $I.$EPrevious error from SILUP$ 401 * occurred in dimension $I.$E'/ 402 data MTXTAB/'LUP(1:$M):$B'/ 403 data MTXTAC/'NDEG(1:$M):$B'/ 404 data MTXTAD/'NTAB(1:$M):$B'/ 405 data MTXTAE/'IOPT(1:$M):$B'/ 406 data MTXTAF/'X((1:$M):$B'/ 407c 408 data MLOC /LTXTAB, LTXTAC, LTXTAD, LTXTAE, LTXTAF, LTXTAG, LTXTAH, 409 1 LTXTAI, LTXTAJ, LTXTAK, LTXTAL, LTXTAM/ 410c 411c 1 2 3 4 5 6 7 8 9 10 412 data MACT / MEEMES,0,0,0, MECONT, MEMDA1,0, METEXT, MEIVEC, 0, 413 1 MECONT, MERET / 414 data INTCHK(0) / 120 / 415 data ID4 / 4 / 416c 417c ***************** Start of executable code ******************* 418c 419 LI = 2 420 NDIMI = NDIM 421 NTABXT = NDIMI+1 422 NTABM = NTABXT + NDIMI 423 LT = NTABM 424 if (NTAB(NTABXT) .le. 0) then 425c Fix NTAB on the first call. 426 if (NDIMI .gt. MAXDIM) then 427 INTCHK(1) = MAXDIM 428 INTCHK(2) = NDIMI 429 IIFLG = -2 430 go to 400 431 end if 432 NTAB(NTABXT) = 1000 433 NTAB(NTABXT+1) = 1 434 do 40 I = 1, NDIMI 435 if (LUP(I) .ge. 4) then 436c LUP(I) is out of range -- fatal error. 437 INTCHK(1) = I 438 INTCHK(2) = LUP(I) 439 IIFLG = -3 440 go to 400 441 end if 442 L = NTAB(I) 443 if (L .gt. 0) then 444c Table is not ragged up to this point. 445 if (I .eq. NDIMI) go to 50 446c If table is ragged save pointer to ragged info. 447 if (NTAB(NDIMI) .lt. 0) NTAB(NTABM+I) = NTABM+NDIMI 448c Get pointer to start of XT data for next dimension 449 if (LUP(I) .eq. 3) then 450 NTAB(NTABXT+I+1) = NTAB(NTABXT+I) + 2 451 else 452 NTAB(NTABXT+I+1) = NTAB(NTABXT+I) + L 453 end if 454 else if ((L .eq. 0) .or. ((L .ne. 1-I) .and. 455 1 ((I .ne. NDIMI) .or. (L .le. -I)))) then 456c Problem in specifying raggedness. 457 LT = NDIMI 458 INTCHK(1) = I 459 IIFLG = -4 460 go to 400 461 else 462 J = NTAB(NTABM+I-1) 463 if (NTAB(I-1) .gt. 0) then 464 NTAB(NTABXT) = -NTAB(I) 465c Get K = number of NTAB entries of extra data 466 K = NTAB(1) 467 do 10 L = 2, -L 468 K = K * NTAB(L) 469 10 continue 470 else 471 K = NTAB(J-1) 472 end if 473 if (NTAB(J) .ne. 0) then 474 LT = J 475 INTCHK(1) = J 476 INTCHK(2) = NTAB(J) 477 IIFLG = -5 478 go to 400 479 end if 480c Change data to be the partial sum of the original data. 481 LT = J + K 482 do 30 L = J + 1, LT 483 if (NTAB(L) .le. 0) then 484c Error, bad value inside ragged table info. 485 INTCHK(1) = L 486 INTCHK(2) = NTAB(L) 487 IIFLG = -6 488 go to 400 489 end if 490 NTAB(L) = NTAB(L) + NTAB(L-1) 491 30 continue 492 if (I .eq. NDIMI) then 493 LT = LT + 1 494 if (NTAB(LT) .eq. -1) go to 50 495c Error -- No end tag where needed. 496 INTCHK(1) = LT 497 INTCHK(2) = NTAB(LT) 498 IIFLG = -7 499 go to 400 500 end if 501c Save index of 0th NTAB entry for extra data for next dim. 502 NTAB(NTABM+I) = J + K + 1 503c Get pointer to start of XT data for next dimension 504 if (LUP(I) .eq. 3) then 505 NTAB(NTABXT+I+1) = NTAB(NTABXT+I) + 2*K 506 else 507 NTAB(NTABXT+I+1) = NTAB(NTABXT+I)+NTAB(J+K) 508 end if 509 end if 510 40 continue 511 end if 512 50 continue 513 IIFLG = 0 514 LINFO(0) = 0 515 LINFO0(0) = 0 516 LINFO0(1) = 1 517 do 60 I = 2, NDIMI 518 LINFO0(I) = LINFO0(I-1) + NDEG(I-1) + 2 519 60 continue 520c Initialize option flags 521 GETERP = .false. 522 SETEXT = 0 523 KAOS = 0 524 IEEOPT = 0 525 STOD = 0 526 LDERIV = 0 527 IOP2N = 254 528 IRAG = NTAB(NTABXT) 529 INTCHK(4) = 0 530 INTCHK(6) = 0 531 INTCHK(5) = -2 * (LINFO0(NDIMI) + 1) 532 LCKEND = 7 533 go to 70 534c 535 62 INTCHK(LCKEND+2) = 1 536 64 INTCHK(LCKEND+1) = IOPT(LI) 537 65 INTCHK(LCKEND) = LOPT 538 LCKEND = LCKEND + 3 539 70 continue 540 LI = LI + 1 541 LOPT = IOPT(LI) 542 if (LOPT .ne. 0) then 543 go to (310, 320, 330, 300, 340, 360), LOPT 544 INTCHK(1) = LI 545 INTCHK(2) = LOPT 546 IIFLG = -8 547 go to 400 548 end if 549 550 INTCHK(1) = LCKEND 551 INTCHK(2) = IOPT(2) 552 INTCHK(3) = 1 553 INTCHK(LCKEND) = 1 554 call OPTCHK(INTCHK, IOPT, 'SILUPM / EOPT$E') 555 if (INTCHK(1) .lt. 0) then 556 IOPT(1) = -12 557 return 558 end if 559 J = LCKEND 560c Store values for storage locations determined by OPTCHK. 561 75 J = J + 1 562 K = INTCHK(J) 563 if (INTCHK(K) .eq. 0) then 564 LERTMP = INTCHK(INTCHK(LCKEND+1)+1) 565 else 566 STOD = INTCHK(INTCHK(LCKEND+2)+1) 567 end if 568 if (J .ne. INTCHK(LCKEND)) go to 75 569 IDXBAS = LERTMP+1 570 IYBAS = IDXBAS + LINFO0(NDIMI) 571 if (STOD .ne. 0) then 572 K = IYBAS + LINFO0(NDIMI) 573 do 80 J = 1, NDIMI 574 KPTD(J) = KPTD(J) + K 575 80 continue 576 LDERIV = KPTD(NDIMI-1) 577 end if 578 IYSAV = IYBAS 579 do 90 I = LERTMP, IYBAS 580c Ensure in SILUP, YT will be defined. (It's ref. but not used.) 581 EOPT(I) = 0 582 90 continue 583 ERRSAV(1) = 0.E0 584 IYT = IDXBAS 585 586 GETERR = GETERP 587 KDIM = 1 588 IXT = 1 589 LUPC = LUP(1) 590 IP(1) = 0 591 100 NT = NTAB(KDIM) 592 110 LEXIT = 0 593c Setup for variable order. 594 if (KAOS .ne. 0) then 595 KAOS = 1 596 IOPTI(3) = NDEG(KDIM) 597 end if 598c Set KEXTRP to indicate how to extrapolate 599 KEXTRP = -1 600 if (SETEXT .ne. 0) KEXTRP = IOPT(SETEXT + KDIM) 601 if (KDIM .ne. NDIMI) go to 180 602c Interpolate in the last dimension. 603c Flag that want results on the exit 604 LEXIT = 1 605c Setup for getting derivatives. 606 if (LDERIV .gt. 0) LEXIT = IOPT(IMAXD+KDIM) + 1 607c Indicate where error info. on Y is found (if any) 608 LEXERR = IEEOPT 609c Indicate whether there are bad points. 610 IOPTI(2) = IOP2N 611c 612 call SILUP(X(KDIM), EOPT(IYSAV), NT, XT(IXT), YT(IP(KDIM)+1), 613 1 NDEG(KDIM), LUPC, IOPTI, EOPT) 614 if (LUPC .lt. 0) LUP(KDIM) = LUPC 615 if (KAOS .lt. 0) go to 390 616 IIFLG = max(IIFLG, IOPTI(1)) 617 130 continue 618c Back up to lower dimension 619 KDIM = KDIM - 1 620 if (KDIM .eq. 0) then 621 Y = EOPT(IYSAV) 622 IOPT(1) = IIFLG 623 if (IOPTI(1) .ge. 0) return 624 IIFLG = IOPTI(1) 625 ERRDAT(1) = EOPT(1) 626 ERRDAT(2) = EBND 627 go to 400 628 end if 629 if (GETERP) ERRSAV(KDIM) = max(ERRSAV(KDIM), EOPT(1)) 630 NDEGEC = LNDEG(KDIM) 631 if (KAOS .ne. 0) then 632 if (LINFO(KDIM) .ge. LINFO0(KDIM) + 3) then 633 NDEGEC = LINFO(KDIM) - LINFO0(KDIM) 634 go to 140 635 end if 636 end if 637 if (LINFO(KDIM) .lt. LINFO0(KDIM) + NDEGEC) go to 220 638c Restore information in the common block 639 140 KGOC = abs(LKGOC(KDIM)) 640 LICINF = LINFO0(KDIM) 641 LINFO(KDIM+1) = LINFO0(KDIM+1) 642 do 150 I = 0, NDEGEC 643 IDX(I) = LIDX(I+LICINF) 644 DX(I) = EOPT(I+IDXBAS+LICINF) 645 150 continue 646 IYSAV = IYBAS + LINFO(KDIM-1) 647 IYT = IYBAS + LINFO0(KDIM) 648c Flag that already have the Y's in order. 649 LUPC = 4 650c Flag that want results on the exit 651 LEXIT = 1 652c Setup for getting derivatives. 653 if (LDERIV .ne. 0) then 654 if (LDERIV .ge. -KDIM) then 655 LDERIV = KPTD(KDIM-1) 656 if (KDIM .ne. 1) LDERIV = LDERIV + NUMD(KDIM) * 657 1 (LINFO(KDIM-1)-LINFO0(KDIM-1)) 658 LEXIT = IOPT(IMAXD+KDIM) + 1 659 end if 660 end if 661c Indicate where error info. on Y is found (if any) 662 if (GETERP) then 663 LEXERR = LERTMP 664 EOPT(LEXERR) = ERRSAV(KDIM) 665 EOPT(LEXERR+1) = 0.E0 666 end if 667c Setup for variable order. 668 if (KAOS .ne. 0) then 669 KAOS = 1 670 IOPTI(3) = NDEGEC 671 end if 672c Indicate no bad points. 673 180 IOPTI(2) = 254 674c 675 if (mod(NDEG(KDIM), 2) .eq. 1) NDEGEC = NDEG(KDIM) 676 call SILUP(X(KDIM), EOPT(IYSAV), NT, XT(IXT), EOPT(IYT), 677 1 NDEG(KDIM), LUPC, IOPTI, EOPT) 678 if (KAOS .lt. 0) go to 390 679 IIFLG = max(IIFLG, IOPTI(1)) 680 if (LEXIT .eq. 0) then 681c Just selected points for this dimension 682 if (LUPC .lt. 0) LUP(KDIM) = LUPC 683c Save info. from the common block. 684 LKGOC(KDIM) = KGOC 685 LNDEG(KDIM) = NDEGEC 686 LICINF = LINFO0(KDIM) 687 LINFO(KDIM) = LICINF 688 do 210 I = 0, NDEGEC 689 LIDX(I+LICINF) = IDX(I) 690 EOPT(I+IDXBAS+LICINF) = DX(I) 691 210 continue 692 go to 230 693 else 694c Just got result for this dimension 695c IOPTI(2) = 0, if variable order and need more points. 696 if (IOPTI(2) .eq. 0) go to 220 697 if (LDERIV .gt. 0) then 698 GETERR = .false. 699 ID = KPTD(KDIM) 700 do 212 I = KDIM+1, NDIMI 701 KC(I) = 0 702 212 continue 703 MXD = MTOTD 704 214 IYDSAV = LDERIV + LEXIT - 1 705 L = KDIM + 1 706 215 if ((KC(L) .eq. IOPT(IMAXD+L)) .or. (MXD .eq. 0)) then 707 MXD = MXD + KC(L) 708 KC(L) = 0 709 L = L + 1 710 if (L .le. NDIMI) go to 215 711 go to 218 712 end if 713 KC(L) = KC(L) + 1 714 MXD = MXD - 1 715 LEXIT = min(MXD, IOPT(IMAXD+KDIM)) + 1 716 LDERIV = IYDSAV + 1 717 do 216 I = 0, NDEGEC 718 EOPT(IYT+I) = EOPT(ID + NUMD(KDIM+1)*I) 719 216 continue 720 ID = ID +1 721 KGOC = abs(KGOC) 722 if (mod(NDEG(KDIM), 2) .eq. 1) NDEGEC = NDEG(KDIM) 723c Interpolate in an outer dimension. 724 call SILUP(X(KDIM), EOPT(IYDSAV), NT, XT(IXT), EOPT(IYT), 725 1 NDEG(KDIM), ID4, IOPTI, EOPT) 726 go to 214 727 218 GETERR = GETERP 728 end if 729 go to 130 730 end if 731 220 LINFO(KDIM) = LINFO(KDIM) + 1 732 if (LDERIV .gt. 0) then 733 if (LINFO(KDIM) .eq. LINFO0(KDIM) + NDEGEC) then 734 if (LKGOC(KDIM) .eq. -2) then 735 LDERIV = -KDIM 736 go to 230 737 end if 738 end if 739 if (KDIM .eq. NDIMI-1) then 740 LDERIV = LDERIV + NUMD(NDIMI) 741 else 742 LDERIV = KPTD(NDIMI-1) 743 end if 744 end if 745 230 IX = LIDX(LINFO(KDIM)) 746c Get here after getting next point in this dimension 747c Compute IP() data for getting XT and YT indices 748 if (KDIM .lt. IRAG) then 749 IP(KDIM+1) = NTAB(KDIM+1) * (IX - 1 + IP(KDIM)) 750 else if (KDIM .eq. IRAG) then 751 K = NTAB(NTABM+KDIM) + IP(KDIM) + IX 752 IP(KDIM+1) = NTAB(K-1) 753 IPRAG = NTAB(K) - NTAB(K-1) 754 else if (NTAB(KDIM) .gt. 0) then 755 IP(KDIM+1) = NTAB(KDIM) * IP(KDIM) + (IX -1) * IPRAG 756 else 757 IP(KDIM+1) = NTAB(NTAB(NTABM+KDIM)+IP(KDIM)+IX-1) 758 end if 759 IYSAV = IYBAS + LINFO(KDIM) 760 KDIM = KDIM + 1 761 ERRSAV(KDIM) = 0.E0 762 if ((LEXIT .eq. 0) .or. (NTAB(KDIM) .lt. 0)) then 763 IXT = NTAB(NTABXT+KDIM) 764 LUPC = LUP(KDIM) 765 if (NTAB(KDIM) .ge. 0) go to 100 766 K = -NTAB(KDIM) 767 K = NTAB(NTABM+K) + IP(K) 768 NT = NTAB(K+IX) - NTAB(K+IX-1) 769 if (LUPC .eq. 3) then 770 IXT = IXT+2*(IP(-NTAB(KDIM))+LIDX(LINFO(-NTAB(KDIM)))-1) 771 else 772 IXT = IXT + IP(1-NTAB(KDIM)) 773 end if 774 go to 110 775 end if 776 if (KDIM .eq. NDIM) then 777 LUPC = LUP(KDIM) 778 go to 110 779 end if 780 LINFO(KDIM) = LINFO0(KDIM) 781 go to 230 782c 783c Process various options 784c Specify error in Y. 785 300 LI = LI + 1 786 IEEOPT = IOPT(LI) 787 INTCHK(LCKEND+2) = 2 788 GETERP = .true. 789 go to 64 790c Return an error estimate 791 310 GETERP = .true. 792 go to 70 793c Setup to set extrapolation degree 794 320 SETEXT = LI 795 LI = LI + NDIMI 796 go to 70 797c Setup for computing extra derivatives 798 330 STOD = IOPT(LI+1) 799 INTCHK(LCKEND+1) = STOD 800 IMAXD = LI + 2 801 MTOTD = IOPT(IMAXD) 802 if (MTOTD .ge. LIDXSZ) then 803 INTCHK(1) = LIDXSZ 804 INTCHK(2) = MTOTD 805 IIFLG = -10 806 go to 400 807 end if 808 KL = IOPT(IMAXD+NDIMI) 809 K = min(MTOTD, NDEG(NDIMI)) 810 if ((KL .lt. 0) .or. (KL .gt. K)) then 811 INTCHK(1) = K 812 INTCHK(2) = KL 813 IIFLG = -11 814 go to 400 815 end if 816 do 332 K = 0, MTOTD 817 LIDX(K+1) = min(K, KL) + 1 818 332 continue 819 NUMD(NDIMI) = KL 820 do 338 J = NDIMI-1, 1, -1 821 L = IOPT(IMAXD+J) 822 K = min(MTOTD, NDEG(J)) 823 if ((L .lt. 0) .or. (L .gt. K)) then 824 INTCHK(1) = K 825 INTCHK(2) = L 826 IIFLG = -11 827 go to 400 828 end if 829 do 334 K = MTOTD, L+1, -1 830 LIDX(K+1) = LIDX(K+1) - LIDX(K-L) 831 334 continue 832 do 336 K = 1, MTOTD 833 LIDX(K+1) = LIDX(K) + LIDX(K+1) 834 336 continue 835 KL = min (KL + L, MTOTD) 836 NUMD(J) = LIDX(KL+1) - 1 837 338 continue 838 INTCHK(LCKEND+2) = NUMD(1) 839 if (STOD .le. 0) then 840 INTCHK(LCKEND+1) = -NUMD(1) 841 INTCHK(LCKEND+2) = LI + 1 842 end if 843 KPTD(1) = 0 844 do 339 J = 2, NDIMI 845 KPTD(J) = KPTD(J-1) + NUMD(J) * (NDEG(J-1) + 2) 846 339 continue 847 INTCHK(5) = INTCHK(5) - KPTD(NDIMI) 848 LI = IMAXD + NDIMI 849 go to 65 850c Setup for variable order. 851 340 LI = LI + 1 852 EBND = EOPT(IOPT(LI)) 853 EBNDR = 0.E0 854 KAOS = 1 855 GETERP = .true. 856 go to 62 857c Setup to indicate bad points 858 360 LI = LI + 1 859 BADPT = EOPT(IOPT(LI)) 860 IOP2N = 255 861 go to 62 862c Error Processing 863 390 IIFLG = IOPTI(1) - 20 864 INTCHK(1) = KDIM 865 MACT(4) = MLOC(12) 866 MACT(2) = -KAOS 867 go to 410 868 400 MACT(4) = MLOC(-IIFLG) 869 MACT(2) = 88 870 if (IIFLG .eq. -1) MACT(2) = 25 871 410 IOPT(1) = IIFLG 872 MACT(3) = -IIFLG 873 call SMESS (MACT, MTXTAA, INTCHK(1), ERRDAT) 874 if (IIFLG .lt. -2) then 875 MACT(7) = NDIMI 876 MACT(10) = NDIMI 877 call MESS (MACT(6), MTXTAB, LUP) 878 call MESS (MACT(6), MTXTAC, NDEG) 879 MACT(9) = MEFVEC 880 call SMESS (MACT(6), MTXTAF, NDEG, X) 881 MACT(9) = MEIVEC 882 MACT(7) = LT 883 MACT(10) = LT 884 call MESS (MACT(6), MTXTAD, NTAB) 885 if (LI .gt. 3) then 886 MACT(7) = LI 887 MACT(10) = LI 888 call MESS (MACT(6), MTXTAE, IOPT) 889 end if 890 end if 891 call MESS (MACT(12), MTXTAA, INTCHK(1)) 892 return 893 end 894