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