1      subroutine SILUPMD (NDIM,X,Y,NTAB,XT,YT,NDEG,LUP,IOPT,EOPT)
2c Copyright (c) 2006, Math a la Carte, Inc.
3c>> 2006-05-03 SILUPMD Krogh Started this debug routine
4c
5c--S replaces "?": ?ILUPMD, ?ILUPM, ?MESS
6c
7c Provide debugging information for Multidimensional Polynomial
8c Interpolation with Look Up, silupm.
9c Design/Code by Fred T. Krogh, Math a la Carte, Inc.
10c
11c *****************     Formal Arguments     ***************************
12c
13c These are exactly as for silupm.  See silupm.f for details.
14c
15c *****************     Variables Referenced     ***********************
16c
17c EOPT   Formal argument, see silupm.f.
18c FDAT   Place to store floating point for error messages.
19c I      Temporary index.
20c IOPT   Formal argument, see silupm.f.
21c IN     Array used to track current index in each dimension.
22c INTCHK Array used for printing integers in error messages.
23c IYT    Pointer to YT information for the current output.
24c J      Index of a 0 in NTAB which is followed by user information
25c        specifying the number of data points as a function of
26c        previously used values from XT.
27c K      Temporary index.
28c KNT    Array used to keep count of number of values to output in
29c        each dimension.
30c LR     Arrray used to find ragged table information in each ragged
31c        dimension.
32c LTXTxx Parameter names of this form were generated by PMESS in making
33c        up text for error messages and are used to locate various parts
34c        of the message.
35c LUP    Formal argument, see silupm.f.
36c LX     Array giving start of XT table for each dimension.
37c MACT   Array giving actions for printing error messages, see MESS.
38c MAXDEG Parameter -- Gives maximum degree of polynomial interpolation.
39c MECONT Parameter telling MESS an error message is to be continued.
40c MEEMES Parameter telling MESS to print an error message.
41c MEIVEC Parameter telling MESS to print an integer vector.
42c MEFVEC Parameter telling MESS to print an floating point vector.
43c MEMDA1 Parameter telling MESS that next item is integer data to be
44c        made available for output.
45c MERET  Parameter telling MESS that this ends the error message.
46c MESS   Subroutine for printing error messages.
47c METDIG Parameter for SMESS to set temporily the digits printed.
48c NDEG   Formal argument, see silupm.f.
49c NDIM   Formal argument, see silupm.f.
50c NDIMI  Internal value for NDIM = number of dimensions.
51c NTAB   Formal argument, see silupm.f.
52c NTABM  = NTABXT + NDIMI -- NTAB(NTABM+I) is 1 + the index of the last
53c        word in NTAB required for ragged table storage through dim. I.
54c NTABXT = NDIMI+1 -- NTAB(NTABXT) = index of first ragged table, =1000
55c        if there is no ragged table.  NTABN(NTABXT+I) = base address
56c        accessing XT information.
57c NUP    Array containing 0, unless an index change in this dimension
58c        triggers a change in a ragged pointer.
59c SMESS  Calls MESS and prints floating data in error messages.
60c X      Formal argument, see silupm.f.
61c XT     Formal argument, see silupm.f.
62c Y      Formal argument, see silupm.f.
63c YT     Formal argument, see silupm.f.
64c
65c *****************     Formal Variable Declarations     ***************
66c
67      integer          NDIM, NTAB(*), NDEG(*), LUP(*), IOPT(*)
68      real             X(NDIM), Y, XT(*), YT(*), EOPT(*)
69c
70c *****************    Parameters and Internal Varialbes ***************
71c
72c       MAXDEG and MAXDIM Should agree with values silupm.f
73      integer MAXDEG, MAXDIM
74      parameter (MAXDEG=15, MAXDIM=10)
75
76      integer I, IIFLG, IN(MAXDIM), INTCHK(4), IYT, J, K, KNT(MAXDIM),
77     1  L, LR(MAXDIM-1), LT, LX(MAXDIM), NDIMI, NTABM, NTABXT,
78     2  NUP(MAXDIM)
79c
80c ************************ Error Message Stuff and Data ****************
81c
82c Parameter defined below are all defined in the error message program
83c MESS and SMESS.
84c
85      integer MEMDA1,MECONT,MERET,MENTXT,METEXT,MEIVEC,MEFVEC,METDIG
86      parameter (MEMDA1=27, MECONT =50, MERET=51, MENTXT=23,
87     1  METEXT=53, MEIVEC=57, MEFVEC=61, METDIG=22)
88c
89      integer MLOC(12)
90      integer MACT(4), MACTDG(3), MACTIV(4), MACTIV1(6),  MACTFV(4)
91      real             FDAT(MAXDIM)
92c
93c ********* Error message text ***************
94c ********* Error message text ***************
95c[Last 2 letters of Param. name]  [Text generating message.]
96cAA Inputs to SILUPM, with NDIM = $I, EOPT dim. of $I.$E
97cAB IOPT($I) = $I, is not a valid option.$E
98cAC IOPT($I) = $I, requests an error estimate.$E
99cAD IOPT($I) = $I, specifies$E
100cAE IOPT($I) = $I, Store derivatives at $I, Max. total derivative $C
101c   of $I,$E
102cAF IOPT($I) = $I, requests Abs. & Rel. errors of: $F $F$E
103cAG IOPT($I) = $I, request accuracy of $F$E
104cAH IOPT($I) = $I, skips data points with YT(?) = $F$E
105cAI LUP($I) = $I; it should be < 4.$E
106cAJ Ragged table badly specified in dimension $I.$E
107cAK Start of ragged table, NTAB($I) = $I; it must be 0.$E
108cAL Middle of ragged table, NTAB($I) = $I; it must be >0.$E
109cAM End of ragged table, NTAB($I) = $I; it must be -1.$E
110c   $
111cAN Extrpolation polynomial degrees:$B
112c   $
113cAO max derivatives per X(I):$B
114c   $
115cAP NDEG():$B
116c   $
117cAQ LUP():$B
118c   $
119cAR NTAB():$B
120c   $
121cAS X(): $B
122c   $
123cAT Dim. 1 XT's = $B
124c   $
125cAU Dim. 1 to $I XT's = $B
126c   $
127cAV Dim. $I XT's = $F + k * %F$E
128c   $
129cAW Dim. $I XT's = $B
130c   $
131cAX YT's = $B
132c   $
133cAY Original ragged size specifications for dimenaion $M:$N$B
134c   $
135cAZ Start of XT specifications for each dimension:$B
136c   $
137cBA Start of specifications for each ragged dimension:$B
138c   $
139cBB Final Ragged Specifications:$B
140c   $
141cBC Dim. 1 to $M Indexes:$B
142c   $
143      integer LTXTAA,LTXTAB,LTXTAC,LTXTAD,LTXTAE,LTXTAF,LTXTAG,LTXTAH,
144     * LTXTAI,LTXTAJ,LTXTAK,LTXTAL,LTXTAM,LTXTAN,LTXTAO,LTXTAP,LTXTAQ,
145     * LTXTAR,LTXTAS,LTXTAT,LTXTAU,LTXTAV,LTXTAW,LTXTAX,LTXTAY,LTXTAZ,
146     * LTXTBA,LTXTBB,LTXTBC
147      parameter (LTXTAA=  1,LTXTAB= 53,LTXTAC= 92,LTXTAD=136,LTXTAE=162,
148     * LTXTAF=232,LTXTAG=286,LTXTAH=325,LTXTAI=375,LTXTAJ=408,
149     * LTXTAK=455,LTXTAL=508,LTXTAM=563,LTXTAN=  1,LTXTAO=  1,
150     * LTXTAP=  1,LTXTAQ=  1,LTXTAR=  1,LTXTAS=  1,LTXTAT=  1,
151     * LTXTAU=  1,LTXTAV=  1,LTXTAW=  1,LTXTAX=  1,LTXTAY=  1,
152     * LTXTAZ=  1,LTXTBA=  1,LTXTBB=  1,LTXTBC=  1)
153      character MTXTAA(3) * (205)
154      character MTXTAB(1) * (34)
155      character MTXTAC(1) * (27)
156      character MTXTAD(1) * (9)
157      character MTXTAE(1) * (8)
158      character MTXTAF(1) * (9)
159      character MTXTAG(1) * (7)
160      character MTXTAH(1) * (16)
161      character MTXTAI(1) * (22)
162      character MTXTAJ(1) * (28)
163      character MTXTAK(1) * (17)
164      character MTXTAL(1) * (9)
165      character MTXTAM(1) * (57)
166      character MTXTAN(1) * (48)
167      character MTXTAO(1) * (52)
168      character MTXTAP(1) * (30)
169      character MTXTAQ(1) * (23)
170      data MTXTAA/'Inputs to SILUPM, with NDIM = $I, EOPT dim. of $I.$EI
171     *OPT($I) = $I, is not a valid option.$EIOPT($I) = $I, requests an e
172     *rror estimate.$EIOPT($I) = $I, specifies$EIOPT($I) = $I, Store der
173     *ivatives at $I, Max.',' total derivative of $I,$EIOPT($I) = $I, re
174     *quests Abs. & Rel. errors of: $F $F$EIOPT($I) = $I, request accura
175     *cy of $F$EIOPT($I) = $I, skips data points with YT(?) = $F$ELUP($I
176     *) = $I; it should be < 4.$ERag','ged table badly specified in dime
177     *nsion $I.$EStart of ragged table, NTAB($I) = $I; it must be 0.$EMi
178     *ddle of ragged table, NTAB($I) = $I; it must be >0.$EEnd of ragged
179     * table, NTAB($I) = $I; it must be -1.$E '/
180      data MTXTAB/'Extrpolation polynomial degrees:$B'/
181      data MTXTAC/'max derivatives per X(I):$B'/
182      data MTXTAD/'NDEG():$B'/
183      data MTXTAE/'LUP():$B'/
184      data MTXTAF/'NTAB():$B'/
185      data MTXTAG/'X(): $B'/
186      data MTXTAH/'Dim. 1 XT''s = $B'/
187      data MTXTAI/'Dim. 1 to $I XT''s = $B'/
188      data MTXTAJ/'Dim. $I XT''s = $F + k * %F$E'/
189      data MTXTAK/'Dim. $I XT''s = $B'/
190      data MTXTAL/'YT''s = $B'/
191      data MTXTAM/'Original ragged size specifications for dimenaion $M:
192     *$N$B'/
193      data MTXTAN/'Start of XT specifications for each dimension:$B'/
194      data MTXTAO/'Start of specifications for each ragged dimension:$B'
195     */
196      data MTXTAP/'Final Ragged Specifications:$B'/
197      data MTXTAQ/'Dim. 1 to $M Indexes:$B'/
198c ********* End of code generated by PMESS **************
199
200c
201      data MLOC /LTXTAA, LTXTAB, LTXTAC, LTXTAD, LTXTAE, LTXTAF, LTXTAG,
202     1  LTXTAH, LTXTAI, LTXTAJ, LTXTAK, LTXTAL /
203      data MACTDG / METDIG, 6, MERET /
204c                      1  2       3      4
205      data MACT / MENTXT, 0, METEXT, MERET /
206c                        1       2  3      4
207      data MACTIV / METEXT, MEIVEC, 0, MERET /
208c                         1  2       3       4  5      5
209      data MACTIV1 / MEMDA1, 0, METEXT, MEIVEC, 0, MERET/
210      data MACTFV / METEXT, MEFVEC, 0, MERET /
211c
212c *****************     Start of executable code     *******************
213c
214      NDIMI = NDIM
215      NTABXT = NDIMI+1
216      NTABM = NTABXT + NDIMI
217      I = 0
218      INTCHK(1) = NDIMI
219      INTCHK(2) = IOPT(2)
220      MACT(2) = MLOC(1)
221      call SMESS(MACTDG, MTXTAA, INTCHK, FDAT)
222 10   continue
223      call SMESS(MACT, MTXTAA, INTCHK, FDAT)
224      if (I .eq. 0) then
225        if ((NDIMI .gt. MAXDIM) .or. (NDIMI .le. 0)) stop
226     1    'Code requires 0 < NDIM < 11.'
227        if (NTAB(NTABXT) .le. 0) then
228c This block of code copied indented, else unchanged from silupm.f,
229c except for code between lines starting with c###
230          NTAB(NTABXT) = 1000
231          NTAB(NTABXT+1) = 1
232          do 40 I = 1, NDIMI
233            if (LUP(I) .ge. 4) then
234c                       LUP(I) is out of range -- fatal error.
235              INTCHK(1) = I
236              INTCHK(2) = LUP(I)
237              IIFLG = -3
238            end if
239            L = NTAB(I)
240            if (L .gt. 0) then
241c                             Table is not ragged up to this point.
242              if (I .eq. NDIMI) go to 50
243c                       If table is ragged save pointer to ragged info.
244              if (NTAB(NDIMI) .lt. 0)  NTAB(NTABM+I) = NTABM+NDIMI
245c       Get pointer to start of XT data for next dimension
246              if (LUP(I) .eq. 3) then
247                NTAB(NTABXT+I+1) = NTAB(NTABXT+I) + 2
248              else
249                NTAB(NTABXT+I+1) = NTAB(NTABXT+I) + L
250              end if
251            else if ((L .eq. 0) .or. ((L .ne. 1-I) .and.
252     1          ((I .ne. NDIMI) .or. (L .le. -I)))) then
253c                                 Problem in specifying raggedness.
254              LT = NDIMI
255              INTCHK(1) = I
256              IIFLG = -4
257              go to 400
258            else
259              J = NTAB(NTABM+I-1)
260              if (NTAB(I-1) .gt. 0) then
261                NTAB(NTABXT) = -NTAB(I)
262c       Get K = number of NTAB entries of extra data
263                K = NTAB(1)
264                do 20 L = 2, -L
265                  K = K * NTAB(L)
266 20             continue
267              else
268                K = NTAB(J-1)
269              end if
270              if (NTAB(J) .ne. 0) then
271                LT = J
272                INTCHK(1) = J
273                INTCHK(2) = NTAB(J)
274                IIFLG = -5
275                go to 400
276              end if
277c       Change data to be the partial sum of the original data.
278              LT = J + K
279
280
281c### This block of code added to that from silupm
282              MACTIV1(2) = I
283              MACTIV1(5) = LT - J + 1
284              call MESS(MACTIV1, MTXTAM, NTAB(J))
285c### End of added block
286
287
288              do 30 L = J + 1, LT
289                if (NTAB(L) .le. 0) then
290c                     Error, bad value inside ragged table info.
291                  INTCHK(1) = L
292                  INTCHK(2) = NTAB(L)
293                  IIFLG = -6
294                  go to 400
295                end if
296                NTAB(L) = NTAB(L) + NTAB(L-1)
297 30           continue
298              if (I .eq. NDIMI) then
299                LT = LT + 1
300                if (NTAB(LT) .eq. -1) go to 50
301c                    Error -- No end tag where needed.
302                INTCHK(1) = LT
303                INTCHK(2) = NTAB(LT)
304                IIFLG = -7
305                go to 400
306              end if
307c       Save index of 0th NTAB entry for extra data for next dim.
308              NTAB(NTABM+I) = J + K + 1
309c       Get pointer to start of XT data for next dimension
310              if (LUP(I) .eq. 3) then
311                NTAB(NTABXT+I+1) = NTAB(NTABXT+I) + 2*K
312              else
313                NTAB(NTABXT+I+1) = NTAB(NTABXT+I)+NTAB(J+K)
314              end if
315            end if
316 40       continue
317        end if
318 50     continue
319c End of code included from SILUPM
320        I = 2
321
322      else if ((INTCHK(2) .le. 0) .or. (INTCHK(2) .gt. 6)) then
323c               A bad option value
324        stop
325      else if (INTCHK(2) .eq. 2) then
326c               Print out desired polynomial degrees for extrpolation.
327        MACTIV(3) = NDIMI
328        call MESS(MACTIV, MTXTAB, IOPT(I+1))
329        I = I + NDIMI
330      else if (INTCHK(2) .eq. 3) then
331c                Print out the derivative information
332        MACTIV(3) = NDIMI
333        call MESS(MACTIV, MTXTAC, IOPT(I+3))
334        I = I + NDIMI + 2
335      end if
336 100  continue
337      I = I + 1
338      INTCHK(1) = I
339      INTCHK(2) = IOPT(I)
340      MACT(2) = MLOC(INTCHK(2) + 2)
341      if (INTCHK(2) .eq. 0) go to 200
342c             1   2    3    4    5    6
343      go to (10, 10, 130, 140, 140, 140), INTCHK(2)
344      MACT(2) = MLOC(2)
345      go to 10
346 130  INTCHK(3) = IOPT(I+1)
347      INTCHK(4) = IOPT(I+2)
348      go to 10
349 140  FDAT(1) = EOPT(IOPT(I+1))
350      if (INTCHK(2) .eq. 4) FDAT(2) = EOPT(IOPT(I+1)+1)
351      I = I + 1
352      go to 10
353
354 200  continue
355      MACTIV(3) = NDIMI
356      call MESS(MACTIV, MTXTAD, NDEG)
357      call MESS(MACTIV, MTXTAE, LUP)
358      MACTIV(3) = NDIMI + 1
359      call MESS(MACTIV, MTXTAF, NTAB)
360      MACTIV(3) = NDIMI
361      if (NTAB(NTABXT) .eq. 0) go to 300
362      call MESS(MACTIV, MTXTAN, NTAB(NTABXT+1))
363      if (NTAB(NTABXT) .lt. 100) then
364        call MESS(MACTIV, MTXTAO, NTAB(NTABXT+NDIMI+1))
365        MACTIV(3) = LT - NTABXT - 2*NDIMI
366        call MESS(MACTIV, MTXTAP, NTAB(NTABXT+2*NDIMI+1))
367      end if
368      MACTFV(3) = NDIMI
369      call SMESS(MACTFV, MTXTAG, INTCHK, X)
370      MACTFV(3) = NTAB(1)
371      call SMESS(MACTFV, MTXTAH, INTCHK, XT)
372c
373c Set up for output of XT and YT
374      do 220 I = 1, NDIMI
375        NUP(I) = 0
376        IN(I) = 1
377        LX(I) = NTAB(NTABXT+I)
378        if (NTAB(I) .lt. 0) then
379          LR(I) = NTAB(2*NDIMI+I)
380          KNT(I) = NTAB(LR(I)+1) - NTAB(LR(I))
381        else
382          KNT(I) = NTAB(I)
383        end if
384        FDAT(I) = XT(LX(I))
385 220  continue
386      IYT = 1
387 240  continue
388c  Print XT, YT data for the current selection.
389      I = NDIMI
390      MACTIV1(2) = I - 1
391      MACTIV1(5) = I - 1
392      call MESS(MACTIV1, MTXTAQ, IN)
393      MACTFV(3) = I - 1
394      INTCHK(1) = I - 1
395      call SMESS(MACTFV, MTXTAI, INTCHK, FDAT)
396      MACTFV(3) = KNT(I)
397      INTCHK(1) = I
398      if (LUP(I) .eq. 3) then
399        call SMESS(MACT(3), MTXTAJ, INTCHK, XT(LX(I)))
400      else
401        call SMESS(MACTFV, MTXTAK, INTCHK, XT(LX(I)))
402      end if
403      call SMESS (MACTFV, MTXTAL, INTCHK, YT(IYT))
404      IYT = IYT + MACTFV(3)
405 250  continue
406      if (NTAB(I) .lt. 0) NUP(-NTAB(I)) = I
407      FDAT(I) = XT(LX(I))
408      I = I - 1
409c               Test if we are done.
410      if (I .eq. 0) go to 300
411      if (NUP(I) .ne. 0) then
412        K = NUP(I)
413        if (LUP(K) .eq. 3) then
414          LX(K) = LX(K) + 2
415        else
416          LX(K) = LX(K) + KNT(K)
417        end if
418        FDAT(K) = XT(LX(K))
419        LR(K) = LR(K) + 1
420        KNT(K) = NTAB(LR(K)+1) - NTAB(LR(K))
421      end if
422      if (IN(I) .lt. KNT(I)) then
423        if (LUP(I) .eq. 3) then
424          FDAT(I) = XT(LX(I)) + IN(I) * XT(LX(I)+1)
425        else
426          FDAT(I) = XT(LX(I) + IN(I))
427        end if
428        IN(I) = IN(I) + 1
429        go to 240
430      end if
431      IN(I) = 1
432      go to 250
433c
434c Restore default for number of digits before exit
435 300  continue
436      call SMESS(MACTDG, MTXTAA, INTCHK, FDAT)
437      print '(''End of output'')'
438      return
439
440
441
442c More error message processing for errors that would show in silupm.f
443 400  continue
444c  Restore state in NTAB.
445      if (NTAB(NTABXT) .lt. 0) then
446        do 410 I = 3*NDIMI+2, LT - 2
447          NTAB(I) = NTAB(I+1) - NTAB(I)
448 410    continue
449      end if
450      NTAB(NTABXT) = 0
451      MACT(2) = MLOC(7-IIFLG)
452      call MESS (MACT, MTXTAA, INTCHK(1))
453      I = 2
454      go to 100
455      end
456