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